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evaluated  usin§  the  experimental  breakdown  parameter^as  proposed  by 
G.  A.  Bird. 'ABeyond  this  limit,  the  molecular  Direct-Simulation-Monte-Carlo 
method  (DSMC)  is  applied..^ 

The  axisymmetric  Monte  Carlo  Simulation  begins  outside  the  nozzle,  where 
the  breakdown  parameter  has  a  value  of  approximately  0.05.  The  actual  shape 
of  this  breakdown  limit  is  a  closed  lobe  surface  which  starts  at  the  nozzle 
lips,  approximately  follows  a  specific  Mach  wave  in  the  Prandtli^t-leyer  fan 
and  closes  back  to  the  axis  far  downstream.  Because  our  interest  is  limited 
to  the  close  vicinity  of  the  corners,  there  the  shape  of  this  limit  may  be 
approximated  to  a  straight  line  (for  an  axisymmetric  flow  making  a  cone). 

For  the  simulation  purposes,  the  domain  between  the  breakdown  limit  and 
the  wall  is  divided  into  sectors,  each  sector  divided  into  radial  regions  and 
each  region  into  simulation  cells.  Each  cell  is  filled  with  a  number  of 
simulated  molecules  relative  to  the  cell  volume  and  local  number  density. 

The  simulation  is  performed  for  each  region  separately  and  contains: 

*  molecular  motions 

*  generation  of  new  molecules  to  simulate  input  flows 

*  deactivation  of  molecules  to  simulate  output  flows 

*  molecular  collisions 

Because  there  is  no  apriori  information  about  the  flow  interaction  between 
different  regions,  the  whole  simulation  is  done  on  an  iterative  basis.  A 
first  run  will  supply  the  output  flux  from  each  region.  These  results  are 
used  as  input  data  for  consecutive  runs. 

'The  program  runs  as  far  as  the  collision  frequency  is  still  high  and 
the  mean  free  path  is  low  compared  with  the  size  of  the  cells.  Beyond 
this  limit  the  flow  may  be  regarded  as  collisionless,  and  the  flux  towards 
the  solid  wall  may  be  computed  directly. 

t 


N  0102-  LF-  OU.  6601 


UnCLASSIFIED _ 

secuniTv  CLASsiriCATiON  or  this  r  AcerwhM  o««  Bmmd) 

i  i 


ACKNOWLEOGtlENT 


Th*  author  wlshea  Co  acknowledge  the  valuable  help  and  cooperation  chat 
Professor  Allen  E.  Fuhs  has  given  him  throughout  this  work. 


DTIC 

Selects 

MAY?  n88 

B 

This  report  has  two  contract  numbers. 

The  one  on  the  back  of  the  front  cover  is 
the  NFS  contract  number.  The  me  on  the 
DD  Form  1473  is  the  DARPA  contract  nuitber. 
Use  both  of  the®  for  the  report. 

Per  Ms.  Rachelle  Paries,  NPS/Code  012 


Accession  For 

NTIF  GR.'.&I 

DTIC  T'.B 

□ 

U  ^inno'inccd 

□ 

j  JtJGtif  1  -  tic.’-- 

- . 

!  By  -  - 

'  Distribution/ _  _ 

1  Avalloh- llty  rod9^ 
n-J-ill  PTii/Or 
iDlat  '  ypooial 


iil 


A.l  DIFFERENT  REGIONS  IN  THE  JET  . 

A. 2  PROGRAM  FLOWCHART  . 

A. 3  PROGRAM  AXSYM  LISTING . 

A. 4  LIST  OF  SYMBOLS  IN  AXSYM  PROGRAM  .  .  . 

A.  5  AXSYM  PROGRAM  -  USER'S  GUIDE  . 

APPENDIX  B:  SIMUL  PROGRAM  . 

B. l  DATA  ORGANIZATION  . 

B.2  MOLECULAR  SIMULATION  FOR  A  GIVEN  REGION 

B.3  SIMUL  PROGRAM  FLOWCHART  . 

B.4  PROGRAM  LISTING  .  . 

B.5  SIMUL  PROGRAM  -  USER'S  GUIDE  . 

B.6  THE  INFLUENCE  OF  THE  Al-IBIENT  GAS  .  .  . 

SUl-lMARY  OF  REPORT . 

LIST  OF  REFERENCES  . 

INITIAL  DISTRIBUTION  LIST  . 


LIST  OF  FIGURES 


1.  THE  SPACECRAFT  LASER  .  2 

2.  THE  HODOGRAPH  PLANE . 9 

3.  FLOW  AT  EXIT  OF  A  SIMPLE  UNDEREXPANDED  JET . 12 

4.  THE  HODOGRAPH  PLANE  FOR  A  CRITICAL  UNDEREXPANDED  JET . 14 

5.  DEPENDENCE  OF  CRITICAL  VALUES  OF  Pamb/Po  0^^  llACH  NUMBER  AT  THE 

EXIT  plane  (Mq)  for  DIFFERENT  VALUES  OF  SPECIFIC  HEAT  RATIO  ( y)  •  •  •  15 

6.  THE  HODOGRAPH  PLANE  FOR  A  HIGHLY  UNDEREXPANOED  JET . 16 

7.  HODOGRAPH  PLANE  SHOWING  THE  POINTS  ON  THE  MESH  Of  CHARACTERISTICS 

(RELATED  TO  PHYSICAL  PLANE  FIGURE  8)  .  19 

8.  SCHEMATIC  SHAPE  OF  A  HIGHLY  UNDEREXPANDED  JET  .  20 

9.  THE  CALCULATION  OF  v  AND  9  FOR  POINT  (3)  IS  BASED  ON  DATA  FOR 

POINTS  1  AND  2 . 21 

10.  CALCULATION  OF  v  AND  0  FOR  AXISYMMETRIC  FLOW . 22 

11.  LIMITS  FOR  CONTINUUM  APPROACH . 25 

12.  REGIONS  IN  A  RING  JET . 29 

13.  CROSS  SECTION  OF  THE  MOLECULAR  SDIULATION  DOMAIN,  DEFINITION  OF 

SECTORS,  REGIONS,  CELLS  AND  COORDINATES  .  33 

14.  VARIATION  OF  THE  ANGLE  UFI . 35 

15.  DEFINITION  OF  INPUT  AND  OUTPUT  FLOWS  OF  SPECIES  1  TO  A  REGION  (KR)  IN 

A  SECTOR  (KS) . 37 

16.  THE  THREE  REGIONS  IN  AN  UNDEREXPANDED  JET . 39 

17.  INDEXING  OF  tlESH  POINTS  IN  THE  AXSYM  PROGKAI-l . 41 

18.  AXSYM  PROGRAM  FLOWCHART  .  42 

19.  ITERATIVE  PROCEDURE  FOR  MACH  NUMBER  CALCULATION  .  44 

20.  THE  MESH  OF  CHARACTERISTICS  IN  REGION  2  AND  RETION  3  -  AXISYMMETRIC 

RING  JET.  Mo-4  ALTITUDE-200km . 45 


2i.  THE  MESH  OF  QIARACTERISTICS  IN  REGION  2  AND  RETION  3  -  AXISYIIMETRIC 


RING  JET.  Mo-2  ALTITUDE-200km . 46 

22.  DEFINITION  OF  SECTOR  GEOMETRY  AND  CELL  VOLUME . 62 

23.  VECTOR  IP(M) . 75 

24.  SIMUL  PROGRAM  FLOWCHART . . . 77 

25.  THE  LOW  DENSITY  REGION  IN  THE  JET . 101 

LIST  OF  TABLES 

1.  ATMOSPHERIC  DATA . . . 6 

2.  CALCULATION  OF  CHARACTERISTICS  ON  PHYSICAL  PLANE  .  18 


ABSTRACT 


Some  procedures  have  been  developed  to  analyze  the  flowfield  of  highly 
underexpanded  axisymmetric  ring  Jets  operated  at  high  altitudes.  The  1-iethod 
of  Characteristics  (HOC)  was  used  to  compute  the  Prandtl-Meyer  expansion  fan 
and  the  flow  parameters  in  that  region*  The  HOC  may  also  be  used  to  obtain 
some  indications  about  the  repetitive  expansion  —  compression  behavior  of  the 
Jet  as  well  as  the  divergent  shape  of  the  expansion  part  downstream,  when 
the  ambient  pressure  goes  below  certain  limits. 

The  continuum  methods  may  be  used  as  far  as  the  limit  where  translational 
equilibrium  ceases  to  exist.  The  breakdown  of  the  continuum  theory  may  be 
evaluated  using  Che  experimental  breakdown  parameter  as  proposed  by 
G.  A.  Bird,  beyond  this  limit,  the  molecular  Dlrect-Simulation-tlonte-Carlo 
method  (OSMC)  is  applied. 

The  axisymmetric  Monte  Carlo  Simulation  begins  outside  the  nozzle,  where 
the  breakdown  parameter  has  a  value  of  approximately  0.05.  The  actual  shape 
of  this  breakdown  limit  is  a  closed  lobe  surface  which  starts  at  tne  nozzle 
lips,  approximately  follows  a  specific  ^iach  wave  in  the  Prandtl-Meyer  fan  and 
closes  back  to  the  axis  far  downstream.  Because  our  interest  is  limited  to 
the  close  vicinity  of  the  corners,  there  the  shape  of  this  surface  may  be 
approximated  by  a  straight  line  (making  a  cone  in  an  axisymmetric  flow). 

For  the  simulation  purposes,  the  domain  between  the  breakdown  surface  and 
the  wall  is  divided  into  sectors,  each  sector  divided  into  radial  regions  and 
each  region  into  simulation  cells.  Each  cell  is  filled  with  a  number  of 
simulated  molecules  relative  to  the  cell  volume  and  local  number  density. 


The  simulation  is  performed  for  each  region  separately  and  contains: 

*  molecular  mo  Cions 

*  generation  of  new  molecules  to  simulate  input  flows 

*  deactivation  of  molecules  to  simulate  output  flows 

*  molecular  collisions 

Because  there  is  no  aprlorl  information  about  the  flow  interaction  between 
different  regions,  the  whole  simulation  is  done  on  an  iterative  basis.  A 
first  run  will  supply  the  output  flux  from  each  region.  These  results  are 
used  as  input  data  for  consecutive  runs. 

The  program  runs  as  far  as  Che  collision  frequency  is  still  high  and  Che 
mean  free  path  is  low  compared  with  the  size  of  the  cells.  Beyond  this  limit 
the  flow  may  be  regarded  as  collisionless,  and  Che  flux  Cowards  Che  solid  wall 
may  be  computed  directly. 


I.  INTRODUCTION 


Gas  jets  released  from  spacecrafts  and  external  flows  about  vehicles  at 
high  altitudes  have  a  renewed  Interest  In  particular  with  regard  to  two  main 
aspects : 

a.  contamination  of  the  spacecraft  walls 

b.  optical  disturbances  caused  by  the  plume. 

A  spacecraft  gas  dynamics  laser  releasing  a  large  quantity  of  gas  Is  highly 
affected  by  these  two  factors  and  the  Interest  In  analyzing  them  and  being 
able  to  control  them  have  a  unique  Importance  In  further  development. 

A  spacecraft  laser  Is  assumed  to  have  a  long  cylindrical  shape  with  the 
optical  power  output  devices  Installed  at  one  of  its  bases  as  shown  in  Figure 
1. 

The  output  gas  is  released  through  a  ring  nozzle,  undergoes  a  fast 
three  dimensional  axisymmetric  expansion,  and  forms  a  plume  covering  the  whole 
meridian  plane  of  the  vehicle.  It  widens  to  large  angles  of  expansion  so  that 
it  may  Intersect  the  laser  beam.  Back  scattered  molecules  may  return  to  the 
wall  of  the  spacecraft  causing  contamination  and  degradation  of  surfaces  and 


vehicle  parts 
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Concinutm  flow  theory  may  be  used  Co  solve  the  flowfield  and  flow 
parameters  as  far  as  there  is  translational  equilibrium,  it  means  that 
intermolecular  interactions  are  fast  enough  to  maintain  expansion  rates. 
Wherever  these  interactions  are  too  slow,  the  continuum  flow  becomes  invalid 
and  the  molecular  flow  theory  should  be  employed. 

The  solution  for  the  continuum  regime  is  computed  here  by  means  of  the 
Method  of  Characteristics  (MOC)  [1,2,3].  The  limit  where  continuum  breakaown 
occurs  was  estimated  by  Che  experimental  breakdown  parameter  as  proposed  by  G. 
A.  Bird  [4].  Beyond  this  limit,  it  is  proposed  to  compute  Che  molecular  flow 
by  means  of  Che  Direct  Simulation  Monte  Carlo  technique  as  described  in  uetail 
by  Bird  [4] . 

For  moderate  and  low  pressure  ratios  (static  pressure  at  the  nozzle  exit 
to  Che  ambient  pressure),  an  underexpanded  jet  exhibits  a  repetitive  expansion 
-  compression  behavior  with  a  geometry  depending  on  the  initial  Mach  angle,  on 
Che  PrandCl-Meyer  fan  angle,  and  on  the  gas  specific  heats  ratio.  For  lower 
ambient  pressure  which  occurs  at  higher  altitudes,  Che  first  compression 
region  is  pushed  out  to  Che  envelope  of  the  jet  forming  the  barrel  shock.  If 
the  ambient  pressure  is  low  enough,  this  compression  region  may  disappear  due 
to  Che  molecular  behavior  of  the  flow  [6,11]. 

The  breakdown  of  the  continuum  theory  occurs  in  a  region  where  the  gas 
density  and  pressure  are  high  compared  with  the  density  and  pressure  in  Che 
ambient  gas.  At  high  altitudes  the  ratios  between  these  parameters  may  reach 
10^  or  more.  Considering  this  range  of  variations,  computational 
validity  dictates  the  use  of  the  Direct  Simulation  Monte  Carlo  method.  In  the 
higher  density  range  the  jet  will  be  considered  as  consisting  of  two  species 
of  gas,  their  molecular  model  will  be  "the  hard  sphere  molecule"  model  and 
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ambient  gas  is  not  allowed  to  protrude.  In  the  lower  density  region  the  flow 
will  be  regarded  as  collisionless. 

In  the  following  chapters  we  bring  the  detailed  description  of  tne 
computer  programs  which  solve  the  different  parts  of  the  flowfleld. 


II.  THE  CONTINUUM  KEGUIE 


A.  THE  TWO  DIMENSIONAL  ISENTKOPIC  UNDEREXPANDED  JET 

The  results  brought  here  are  based  on  the  supersonic  steady  isentroplc 
flow  theory  as  described  in  literature  (see  for  example,  Shapiro  [1],  Llepman 
and  Roshko  [2],  and  Owczarek  [3]). 

The  ranges  of  paremeters  of  a  Jet  emerging  from  a  gas  dynamics 
spacecraft  laser  are: 

a.  The  Mach  number  at  the  exit  surface  Mq  ■  4.  The  static  pressure 
at  the  exit  plane  Pq  ■  136pa.  Ambient  pressure  (Pamb)*  temperature 
(Tamb)  ^nd  other  thermodynamic  properties  of  ambient  gas  depend  on 
the  altitude  as  shown  in  Table  1.  The  Jet  gas  may  consist  of  DF,  HF, 
Helium  and  other  species.  In  the  programs  we  limit  the  composlton  to 
two  species:  Air  and  He.  (The  program  allows  changes  in  the 
composition  and  types  of  gas.) 

P 

b.  The  pressure  ratio  —  for  the  minimum  required  altitude  (20(j 

‘^amb 

km)  assures  that  the  Jet  is  highly  underexpanded  (we  show  later  the 
influence  of  this  ratio  on  the  shape  of  the  Jet). 

The  following  thermodynamic  relations  are  valid  as  long  as  the 
compesslble  flow  is  isentroplc 
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where  Px>  and  px  are  the  total  temperature,  pressure  and  density 

(constant  for  Isentroplc  field).  T,  P,  and  p  are  local  temperature,  pressure 
and  density  y  is  the  specific  heat  ratio  of  the  gas  (considered  here  as 
constant),  ti  is  the  local  Mach  number. 

The  partial  differential  equation  of  motion  for  supersonic  2-D 
Irro rational  and  Isentroplc  flow  is  a  hyperbolic  equation  having  solutions 
obtained  from  invariants  along  characteristic  lines.  Physical  interpretation 
of  these  lines  are  the  compression  or  expansion  waves  which  are  oriented  at 
Mach  angles  relative  to  the  streamlines. 

Unce  the  directions  of  the  characteristics  (waves)  are  determined 
everywhere  in  the  field,  all  other  parameters  may  be  calculated. 

B.  THE  TWO  UUIENSIONAL  PLANAR  JET 

Tne  compressible  supersonic  Jet  flow  is  characterized  by  two  families 
of  characteristics  (pressure  waves)  starting  at  each  corner  of  the  nozzle 
lips.  Each  of  these  families  of  waves  forms  a  Frandtl-lieyer  fan.  The 
streamlines  crossing  the  characteristic  waves  bend  outwards  resulting  in  an 
increase  in  the  flow  area.  The  angle  y  between  the  streamline  and  the 
pressure  wave  is  a  function  of  the  local  Mach  number  as 

y  »  arcsin  (1/M)  (4) 

The  symbol  y  is  the  tlach  angle. 


Using  the  isentropic  relations,  we  can  find  the  relation  between  the 


turning  £mgle  (9)  and  the  local  t-iach  number  (li)  as 


d9  -  - 


(M  1) 


Integration  between  the  conditions  M*1  and  M  gives  the  total  turning 
angle  starting  at  the  throat  (where  M»l)  up  to  a  point  with  given  M.  The 
result  gives  the  Prandtl-Meyer  function  (angle)  as 

v(M)  ■  arctg  (M^-1)  -  arctg  (fI^-1)  (6) 

In  the  close  vicinity  of  the  nozzle  lips  wnere  the  two  families  of 
characteristics  do  not  Intersect  with  each  other  there  is  a  "simple  region"  of 
expansion.  There  the  flow  parameters  are  defined  by  v  and  9  of  each 
characteristic  line.  Further  downstream  the  waves  Intersect  each  other.  In 
this  part  of  the  flow  parameters  are  defined  by  the  two  intersecting 
characteristics.  A  singularity  occurs  when  the  initial  Mach  number  is  unity 
and  a  special  treatment  is  required  to  start  the  calculations  at  that  point 
(this  special  treatment  has  not  been  brought  here).  A  particular  Importance 
of  the  Prandtl-Meyer  function  is  when  analyzing  the  two  dimensional  flow  using 
the  hodograph  plane. 


C.  THE  HODOGRAPH  PLANE  FOR  A  TWO  DIMENSIONAL  JET 

The  hodograph  plane  is  a  representation  of  the  flow  parameters  in  the 
velocity  plane.  Figure  2  shows  a  hodograph  plane  calculatea  for  a  simple 
gas  (yi'l.A).  The  circles  represent  constant  Mach  numbers  and  constant 


The  epicycloids  are  the  two  families  of 


pressure  ratios  (~) 

T 


characteristics.  The  angles  6  are  the  turning  angles  due  to  expansion  or 
compression  (which  are  both  present  in  supersonic  jets). 

To  define  an  isentroplc  supersonic  jet  on  the  hodograph  plane  it  is 
necessary  to  define  the  Mach  number  at  the  exit  plane,  the  pressure  ratios 

p 

(•^)  and  (-5 — -)  ,  and  the  specific  heat  ratio  (y)  of  the  gas. 
i't 

Investigating  the  shapes  of  the  jet  as  a  function  of  ranges  of  parameters 
we  may  get: 


a. 


simple  underexpanded  jets  for  which 


^amb 


is  high  enough  so  that 


Che  two  families  of  characteristics  intersect  each  other. 


b.  a  critical  underexpanded  jet  for  which 


is  low  enough  so  that 


the  Intersection  between  the  outer  characteristic  lines  of  the  two 
familes  occur  on  Che  outer  hodograph  circle. 

c.  highly  underexpanded  jets  for  which  a  part  of  characteristics  do  not 
intersect  at  all. 


d.  expansion  into  complete  vacuum  so  that  there  are  no  reflections  from 
Che  jet  boundaries  and  therefore  the  compression  region  disappears. 


D.  THE  SHAPES  OF  TWO  DIMENSIONAL  JETS 

The  following  paragraphs  further  detail  the  different  shapes. 

1 .  Simple  Underexpanded  Jets 

Figure  3  shows  the  physical  plane  and  the  hodograph  plane  of  a  simple 

underexpanded  jet.  If  Mo>l  ,  an  initial  turn  of  the  flow  0q  Is  made  within 
the  nozzle.  An  aditional  turn  of  6  is  due  to  the  underexpansion.  6  is  found 
by  the  intersection  of  characteristics  (1-2)  with  the  circle  defined  by 

^amb/^T  • 

When  Mo>l  >  the  characteristic  line  1-2  is  described  in  the  physical 
plane  by  a  region  in  which  only  one  family  of  expansion  waves  are  present 
(simple  region,  see  transverse  line  between  1  to  2  in  the  physical  plane).  A 
different  family  of  expansion  waves  forms  a  second  simple  region  when  moving 
between  2  to  3  .  At  a  larger  distance  from  the  exit  plane,  reflected  waves 
from  the  free  streamline  (jet  boundary)  cause  compression.  An  ideal 
representation  of  such  a  jet  is  a  repetitive  pattern  of  expansion  and 
compression. 

For  Mg^l  I  the  tail  waves  of  both  families  are  perpendicular  to 
the  flow,  thus,  both  lie  on  line  AB  .  That  means  that  if  Mq*!  there  is  no 
simple  region  near  the  exit  plane  of  the  jet.  The  characteristic  line  1-2  on 
the  hodograph  plane  becomes  a  single  point  A  (or  B)  when  located  In  the 
physical  plane. 

2 .  Critical  Underexpanded  Jets 

We  define  a  critical  underexpanded  jet  when  point  (3)  (see  Figure  4) 
lies  on  the  limiting  circle  !(„  or  P/F„>0.  This  means  that  there  is  a  core 
within  the  jet  where  the  pressure  approaches  zero  and  Ilach  approaches 
infinity.  This  core  is  tneoretlcally  bounded  at  its  upstream  side  by 


expansion  waves  and  downstream  by  compression  waves 


Figure  3 
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(b)  Hodograph  Plane 


Using  the  theoretical  expressions  for  isentropic  ideal  flow  one 
P  P 

may  derive  the  values  of  — -  or  -  as  function  of  M  which  causes  a  jet 

to  be  critical  underexpanded. 

P 

Figure  5  shows  results  of  p -  as  functions  of  Mg  for  gases  with 

o 

different  specific  heat  ratio. 


3.  Highly  Underexpanded  Jets 

When  Pamh  is  lower  than  the  critical  values  as  shown  in  figure  4, 
point  3  does  not  exist  (there  is  no  Intersection  between  lines  2-3  2 '-3'). 

This  means  that  the  repetitive  reversible  expans  ion/ compress  ion  shapes  ceases 
to  exist.  The  envelope  of  the  Jet  starts  at  an  angle  defined  by  9  at  the 
nozzle  exit  plane  and  approaches  an  asymptotic  angle  defined  by  (see 

Figure  6). 

In  this  case  we  get  two  (symmetric)  groups  of  characteristics  C1-C2  and 
C'|-C'2  (Figure  6)  with  no  intersection  between  them.  and  '  define  the 
inner  limit  for  reflected  characteristics,  C2  and  C2'  define  the  outer  limit. 
As  the  rest  of  reflected  waves  lie  between  C2  or  C'2  and  the  jet  bounaary,  we 
may  conclude  that  a  compression  region  may  exist  only  in  a  layer  along  the  jet 
boundary. 

because  of  irreversible  effects  such  as  shear  stresses,  heat  transfer  due 
to  high  temperature  gradients,  or  condensation  effects  in  real  gases,  the 
compression  layer  may  be  Interpreted  as  the  "barrel  shock". 


.2> 


A 


.5 


J _ L 

J  .2. 


Figure  5.  Dependence  of  cr Ideal  values  of  ^amb^^o 

on  the  Mach  number  at  Che  exit  plane  (M^)  for  different  values 
of  specific  heat  ratio  (y). 
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Figure  6.  The  Hodograph  Plane  for  a  Highly  Underexpanded  Jet 


(y  -  1.4) 

**  naxlmum  turning  angle 

9  **  total  turning  angle  In  the  specific  Jet  6  >  9 

Slim  ~  limiting  angle  of  compression  region. 


Figure  (7)  shows  the  hodograph  plane  for  an  air  jet  (-pi. 4)  expanding 
into  an  ambient  pressure  105  times  lower  than  the  static  pressure  at  the  exit 
plane.  Figure  (8)  is  a  schematic  description  of  the  jet  (the  data  for 
this  jet  is  given  in  Table  2).  The  shape  shown  in  Figure  (8)  may  be  compared 
with  the  barrel  shock  photograph  in  page  208  Reference  [3]. 

4.  EXPANSION  INTO  A  COMPLETE  VACUUM  -  Pamb~Q 

This  is  an  extreme  situation  in  which  the  maximum  turn  angle  of 
streamlines  occur  near  the  nozzle  exit  plane.  The  theoretical  free  stream  is 
defined  only  by  the  theoretical  turning  angles.  All  streamlines  in  the  flow 
expanded  monotonlcally  towards  P^O  without  being  reflected  by  the  jet 


boundaries . 


Figure  8.  Schematic  shape  of  a  highly  underexpanded  jet 
(See  Figure  7  -  The  Hodograph  Plane) 

The  arrows  Indicate  flow  direction. 


E.  THE  HETHOD  OF  CHARACTERISTICS;  COMPUTATION  OF  PLANAR  AND  AEISTMMETRIC 
THO-ODfENSIONAL  FLOWS 

For  a  planar  two  dlnanalonal  flow,  the  PrandCl**Meyer  function  v  and 
flow  direction  9  at  any  point  (3)  in  the  field  nay  be  calculated  using  data  of 
two  other  points  (1)  and  (2)  located  on  characteristic  lines  that  intersect  at 
(3)  (see  Fig.  (9)) 

'*3  "  T  '*‘T 

®3  “  T  “  '*2^  T  ^  ®2^ 


Figure  9.  The  calculation  of  v  and  9  for  point  (3)  is  based  on  data  for 
points  1  and  2. 


2 


For  aidsyiflnetrlc »  two^dlnonsloiuil  flow^  Llepnann  and  Roshko  [21  (teveloped 
expressions  for  finding  the  Prandtl-*Meysr  function  (v)  and  the  flow  direction 
(0)  which  are  given  by: 

1  1  1 

”3  ■  7  <''1  *  '’2>  *  T  '*1  ■  *  7  "l  ‘^13  *  "2  ‘"23* 

(9) 

1  1  I  ®^“®2 

®3  "  I  ■  ''2^  '*‘1  ^®1  ®2^  2  “l  "77”  ^^13  "  “2  "7^  ^’’23^ 

(10) 


The  angles  and  subscripts  are  shown  in  Figure  (10). 


Figure  10.  Calculation  of  6  and  v  for  axisymmetric  flow. 


It  is  obvious  that  for  the  axially  sycmetrlc  flows  the  Increase  in  the 
radius  causes  the  Increase  in  the  flow  cross  section  and  Influences  the  flow 
direction  and  the  Prandtl-Ueyer  function.  These  facts  have  been  taken  into 
account  when  developing  the  "compatibility  equations"  (9,10). 

Using  equations  (9,10)  we  have  developed  a  computer  program  which 
enables  the  calculation  of  the  jet  flow  for  a  two  dimensional  and  for  an 
axially  symmetric  geometry  (ring  Jet). 

The  listing  of  the  program,  the  program  description  and  some  results  are 


given  in  Appendix  A 


III.  THE  BREAKDOWN  OF  THE  CONTINUUM  THEORY 


A.  GENERAL  CRITERIA 

As  described  in  detail  by  Bird  (Chapter  1  Reference  [4]),  the  validity 
of  the  continuum  approach  has  been  identified  with  the  validity  of  the  Navler 
Stokes  equations.  This  requires  that  the  Knudsen  number  Kn  ■  VL  should  be 
small  compared  with  unity  (X  is  the  mean  free  path  and  L  is  a  scale  length  for 
the  specific  flow  field).  For  K^  larger  than  a  certain  limit  (between  0.01 
to  0.1  depending  on  the  required  accuracy)  a  microscopic  approach  is 
necessary. 

For  small  values  of  L  ,  the  microscopic  approach  may  lead  to  statistical 
fluctuations  of  Che  results  due  Co  Che  small  number  of  molecules  participating 
in  the  flow  processes.  In  figure  (11)  which  was  reproduced  from  Bird's  book 
[4,  figure  1.6],  the  regimes  of  rarefied  flow  and  high  fluctuations  are 
depicted.  The  flow  around  Che  Jet-air  boundaries  near  a  spacecraft  is 
generally  rarefied  (high  Knudsen  number),  but  has  insignificant  fluctuations. 


nuirucicrbiic  ilinicniiiiMi  /.  (iiwlroi) 


B.  THE  EMPIRICAL  CRITERION 


The  Method  of  Characteristics  (MOC)  was  used  to  compute  the  jet  flow 
and  the  results  obtained  from  the  computer  program  ”AXSYM"  are  valid  as  long 
as  the  continuum  flow  theory  is  valid. 

The  continuum  flow  requires  that  the  mean  free  path  should  be  negligibly 
small  in  comparison  with  the  scale  length  of  the  macroscopic  flow  variations. 
The  classical  theory  for  Prandtl-Meyer  expansion  may  therefore  be  expected  to 
fail  at  progressively  larger  distances  from  the  nozzle  lips  as  the  gas  density 
decreases  with  the  increasing  flow  angle  and  Mach  number.  The  empirical 
criterion  for  the  breakdown  of  continuum  flow  in  steady  expansion  flow  [4]  is 
that 

P  =  -3-  ^  =  0.05  Ul) 

p  V  dS 

where  q  ■  stream  velocity 

p  »  dens i ty 

V  *  molecular  collision  frequency 

>  absolute  change  in  density  while  moving  a  distance 
I 

dS  along  a  streamline 
Introducing  the  breakdown  parameter  P  into  the  program,  gives  the 
aefinltion  of  the  boundary  where  the  flow  should  be  calculated  by  means  of  the 
molecular  flow  theory,  i.e.,  by  solving  the  Boltzmann  equation. 
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For  an  under expanded  jet  with  a  high  initial  Mach  number, 
the  breakdown  surface  is  nearly  a  streamline.  Furthermore,  the  range  of  flow 
parameters  for  the  present  problem  are  such  that  the  simple  region  extends  to 
very  large  distances  and  near  the  nozzle  lip  the  breakdown  limit  may  be 
approximated  by  a  straight  line. 

For  the  axisymmetric  jet  there  is  no  simple  region,  however,  for  the 
region  of  interest  it  may  be  regarded  as  linear. 

The  method  proposed  in  Che  present  work  for  solving  the  flow  behind  the 
breakdown  boundary  is  Che  Direct  Simulation  Monte  Carlo  (DSMC).  For  tnis 
purpose,  a  computer  program  “SIMUL”  was  developed.  In  the  following  chapter 
we  describe  the  algorithms  required  for  the  specific  problem,  the  geometry  and 
Che  data  organization.  Detailed  program  description  is  given  in  Appendix  (B). 
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THE  MOLECULAR  FLOW  IN  AN  AXISYM1ETRIC  RING  JET 


I  A.  GENERAL  CONSIDERATIONS 

The  part  of  the  field  in  which  the  Jet  may  be  calculated  by  means  of  the 
continuum  theory  was  described  in  Chapter  II.  There  we  calculated  also  the 
I  boundaries  where  continuum  theory  becomes  invalid  and  molecular  calculation 

should  be  employed.  In  fact,  the  molecular  theory  and  the  molecular  Boltzmann 
equations  are  universal  and  hold  for  the  entire  flowfield.  However, 

I  computational  requirements  make  the  Boltzmann  equation  impractical  for  the 

upstream  flows.  Therefore  we  limit  our  solution  only  to  the  part  of  the  flow 
beyond  the  region  where  continuum  breakdown  occurs. 

!  As  a  result  obtained  from  MOC  solution  the  "breakdown”,  i.e.,  the  locus 

where  the  breakdown  parameter  p  has  values  between  0.03  to  0.06,  for  the 
region  close  to  the  nozzle  lips  this  boundary  may  be  approximated  by  a 
I  straight  line  (for  axisymmetrlc  flow  this  line  is  the  envelope  of  a  cone,  see 

Figure  (12)). 

For  the  specific  jet  and  gas,  the  breakdown  occurs  in  a  region  where  the 
^  number  density  is  in  a  range  of  10^1  molecules/m^ .  For  ambient  gas  at  an 

altitude  of  200  km  the  number  density  is  10 and  decreases  to  a  range  of  10 ^ ^ 
at  1000  km.  In  order  to  be  able  to  express  this  vast  change  in  a  simulation, 

i 

'  we  would  need  to  have  an  extremely  large  number  of  cells  which  would  be 

impossible  to  store  in  a  computer.  To  overcome  this  problem  we  are  required 
!  to  make  a  less  exact  formulation  which  enables  the  production  of  results, 

j 


having  to  pay  the  penalty  of  "smearing"  the  steep  gradients  and  obtaining 
averages  within  layers  of  simulated  cells.  Unfortunately  it  is  impossible  to 


predict  how  far  the  simulated  results  will  be  from  the  exact  solutions.  T 
comparisons  have  to  be  made  after  getting  final  results  of  this  simulation 


1 ,  The  Direct  Simulation  Monte  Carlo  Method 

The  direct  simulation  Honte  Carlo  Method  Is  a  technique  for  a  computer 
modeling  of  a  real  gas  by  some  thousands  of  simulated  molecules.  The  velocity 
components  and  the  coordinates  of  the  simulated  molecules  are  scored  In  the 
computer  and  are  modified  with  time  as  a  result  of  collisions  and  boundary 
interactions.  A  detailed  description  of  some  problems  and  their  solutions 
by  means  of  direct  simulation  is  given  in  [4j. 

To  follow  the  molecular  motion  It  is  necessary  to  divide  the  simulated 
domain  into  a  network  of  cells .  The  size  of  ^  cell  must  be  such  that  the 
change  in  flow  properties  across  each  cell  Is  small .  The  time  is  advanced  in 
discrete  steps  DTM,  such  that  DTM  Is  small  compared  with  the  mean  collis Ion 
time  per  molecule.  If  there  is  a  flow  going  through  the  domain,  DTM  should  oe 
small  compared  with  Che  mean  time  required  for  the  mean  flow  to  cross  tne 
cells.*  Both  cell  size  ( ( DR) , ( DDALFA)  ■-  radial  size  and  angular  size  as  they 
appear  In  Che  program)  and  Dill  may  vary  in  Che  simulation  with  position  and 
time. 

Applications  such  as  free  jet  expansion  in  which  large  gradients  of  flow 
properties  are  expected,  may  require  a  very  large  nuraber  of  cells  for  the 
simulation.  In  these  cases  the  computer  memory  requirements  to  store  cells' 
data  and  molecules'  data  may  exceed  the  available  computer  storage.  A 


*If  DTH  is  chosen  to  be  very  small  compared  with  the  mean  time  between 
collisions  then  the  simulation  will  require  a  very  large  number  of  runs  such 
that  the  number  of  collisions  will  be  sufficient.  If  DTM  is  large  the 
molecules  are  washed  out  by  the  mean  flux  and  there  is  no  time  for  the 
collisions  to  influence  tne  flow. 


solution  for  this  problem  is  to  divide  the  simulation  space  into  smaller 
regions  and  to  run  the  simulation  for  each  region  separately.  If  there  is  an 
Interaction  between  different  regions,  which  apriori  is  undefined,  the 
solution  should  be  found  Iteratively.  (That  means  that  each  run  will  provide 
data  for  consecutive  runs  and  the  procedure  should  be  repeated  until  the 
results  converge  to  a  steady  solution. 

The  computation  of  a  representative  set  of  collisions  based  on  mean 
collision  time  per  molecule  is  invalid  for  a  computerized  simulation  because 
of  the  large  computer  time  and  computer  memory  requirements.  Instead,  the 
method  proposed  by  Derzko  which  is  described  in  details  by  Bird  [4]  may  be 
employed.  Following  this  method,  an  averaged  mean  time  between  collisions  of 
species  L  with  species  M  for  a  cell  is  calculated.  The  number  of 
collisions  of  each  type  (L.M  species)  is  such  that  the  collision  time 
counters  are  kept  concurrent  with  the  overall  time  parameter.  The  L.M 
collision  time  for  a  cell  containing  Nl  and  Nji  molecules  with  collision 
cross  section  olh  ,  number  densities  n^  and  nf.i  and  relative  velocity  Cj., 
is  given  by 


At 


c 


M.  1  ^  Mp  1 


(12) 


where  LF  and  HP  are  the  probabilities  that  the  collision  will  be  effective 
for  the  L  and  M  molecules  respectively. 


B.  THE  GEOMETRY  OF  THE  SIMULATED  DOMAIN,  SECTORS,  REGIONS  AND  CELLS 

Figure  (13)  shows  a  cross  section  of  the  simulated  domain  for  the 
axlsymmetrlc  (ring)  flow.  Points  A  and  A'  are  the  nozzle  lips.  Starting  at 
"A"  and  assuming  the  "breakdown"  boundary  to  be  a  straight  line,  we  obtain  the 
cross  section  of  the  molecular  domain  as  a  sector  defined  by  LAM.  The  solid 
wall  is  defined  by  AL.  The  arc  LM  may  be  assumed  to  be  far  enough  so  that  the 
pressure  along  it  may  be  assumed  to  equal  the  ambient  pressure.  Molecules 
originated  in  the  jet  cross  the  breakdown  boundary  with  a  velocity,  direction, 
temperature  (and  other  thermodynamic  properties)  as  found  from  the  continuum 
solution. 

The  molecular  domain  LAM  is  divided  into  secondary  sectors,  and  each  of 
these  are  divided  into  several  radial  regions  making  the  "simulation 
regions " . 

Because  we  have  no  apriori  information  on  how  the  expansion  occurs,  the 
angle  of  each  sector  (which  mainly  is  in  the  direction  of  the  expansion 
gradients)  is  left  to  be  a  result  of  the  internal  calculation. 

Each  region  is  dlvlaed  into  NRD  radial  divisions  and  NAD  angular 
divisions  making  a  network  of  NAD*NRD  simulation  cells.  The  angle  DALFA  of 
all  regions  in  a  sector  is  constant.  Taking  NAD  constant  for  all  regions  in 
a  sector,  we  get  the  angle  of  a  cell  DDALFA  constant.  Defining  the  radial 
size  of  a  cell  DR  as  constant  we  get  a  cell  cross  section  area  proportional 
to  the  radius  R  measured  from  the  nozzle  lip  (point  A). 


The  size  of  a  cell:  In  order  to  get  accurate  simulated  results  it  is 
recommended  to  define  the  size  of  a  cell  (DR  and  R*DOALFA)  small  compared  with 
the  mean  free  path  of  the  molecules  X  (typical  DR»X/3).  However,  as  we  do 
not  expect  to  get  large  changes  in  flow  parameters  along  the  radius  we  may 
allow  DR  be  much  larger  than  X/3.  The  angular  size  of  the  largest  cell  in  a 
sector  should  comply  with  this  requirement,  but  because  of  the  computer 
limitation  it  is  set  to  be  equal  to  5*X.  This  will  be  the  basis  for  aeflning 
DALFA  for  each  sector. 

C.  INITIAL  NUMBER  OF  MOLECULES  IN  CELLS 

A  "reasonable"  number  of  molecules  in  a  simulation  is  several  thousands 
(a  larger  number,  which  is  better,  may  be  used  for  simple  problems  or  when 
using  a  single  user  computer  with  large  user  memory  space).  The  initial 
setting  of  molecules  in  cells  is  usually  based  on  a  guess  of  the  number 
density  in  the  specific  cell.  (The  number  of  molecules  in  cells  will  change 
during  the  simulation  according  to  the  input/output  calculated  fluxes  to  the 
specific  region). 

The  number  density  and  the  size  of  a  cell  are  specified  only  in  three 
dimensional  flows.  When  applied  to  a  two  dimensional  flow  the  simulation  may 
be  regarded  as  applying  to  an  arbitrary  thin  slice  of  the  real  flow.  In  the 
axisymmetric  flow  we  define  the  width  of  a  cell  by  the  angle  DFI  as  shown  in 
Figure  (14),  constant  within  a  region. 

The  initial  number  aensity  in  cells  of  a  given  region  is  set  constant. 
Defining  the  total  number  of  simulated  molecules  in  the  region  the  number  of 
molecules  in  each  specific  cell  becomes  a  function  of  DFI.  For  example, 
assume  that  we  limit  the  number  of  molecules  in  the  smallest  cell  in  the 


region  to  15  then; 


Variation  of  the  angle  DPI. 

To  nalntaln  the  number  of  simulated  molecules  within  computational 
limits,  the  'width*  of  each  region  defined  by  DPI  la  such  that: 

MIN  >  number  of  molecules  In  smallest  cell  In  a  region 
MIN  •  VOLUME  (smallest  cell)  *  number  density 

■  f(R,  ALPA,  DALPA)  *  DPI  *  number  density  for  flux 
calculations 

DPI  la  a  weighting  factor. 


DFI*(contant)  *  (number  density)  ■  15 


Other  cells  contain  the  initial  number  of  molecules  proportional  to  their 
volumes . 

D.  DEFINITION  OF  INPUT  AND  OUTPUT  FLOWS  FOR  A  REGION 

The  cross  sections  of  all  regions  (except  those  regions  near  the  nozzle 
lips)  are  quadrilateral.  Through  the  sides  of  the  region  molecules  are 
allowed  to  enter  or  to  leave  according  with  the  boundary  conditions  or  as  a 
result  of  molecular  velocity.  For  the  first  sector,  near  the  breakdown 
boundary  the  input  flow  (FWPl  and  FWP2  see  Figure  (15))  is  defined  by  the 
results  from  the  continuum  flow.  FENl,  FNNl,  FSNl,  etc.  are  results  of 
counting  and  averaging  the  outgoing  molecules  (the  different  vector  names  will 
be  explained  in  Appendix  B) . 

For  the  neighbor  regions  these  output  fluxes  become  inputs  and  have  to 
be  adjusted  according  to  the  differences  in  the  angle  DFl  of  the  different 
regions . 

The  simulation  starts  with  regions  in  the  sector  near  the 
breakdown  boundary.  At  this  time  there  is  no  data  for  input  flows  through 
faces  E  and  N  of  the  cell.  An  additional  run  of  the  whole  program  is  required 
in  order  to  take  these  calculated  flows  into  account.  If  the  accuracy  of  the 
results  is  Important  we  may  run  this  type  of  iteration  several  times  until 
the  results  become  stable.  (Only  after  running  the  program  tor  the  whole 
domain  once  we  shall  be  able  to  evaluate  the  importance  of  these  iterations.) 
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E.  COLLISIONLESS  FLOW 


In  several  sectors  near  the  breakdown  boundary  we  may  find  a  high  number 
density  and  the  mean  free  path  small  compared  with  the  size  of  a  cell.  There 
the  calculated  collisions  are  expected  to  have  an  Influence  on  the  flow 
parameters.  For  wider  expansion  angles  the  collisions  become  rare  mainly 
because  of  the  decrease  In  the  density.  In  the  ambient  gas  the  mean  free  path 
(for  200  km  altitude)  Is  240  m.  Comparing  this  number  with  the  size  of  the 
simulated  domain  may  lead  to  the  conclusion  that  there  the  flow  may  be 
regarded  as  collisionless. 

We  may  define  a  limiting  line  In  the  flow  where  the  collisions  become 
Insignificant.  Consequently,  molecules  crossing  this  limit  will  In  fact 
continue  moving  In  straight  lines;  a  part  of  them  reach  the  solid  wall. 

Introducing  this  Idea  of  the  collisionless  flow  we  may  reduce  the 
computation  time  and  the  memory  requirements. 

F.  TWO  DIMENSIONAL  PLANAR  FLOW  VS.  AXISYMMETRIC  FLOW 

The  cell  dimensions  are  completely  specified  only  In  three  dimensional 
flow.  When  applied  to  two  dimensional  flow,  the  simulation  may  be  regarded  as 
applying  to  an  arbitrarily  thin  slice  of  the  real  flow.  The  thickness  of  the 
slice  may  be  chosen  such  that  the  number  of  simulated  molecules  complies  with 
the  cell  volume  and  the  physical  number  density.  For  the  axlsymmetrlc  flow  we 
have  defined  the  angle  DFI  as  the  third  coordinate  so  that  the  volume  of  the 
cell  Is  completely  specified. 

Once  the  geometry  Is  defined,  the  simulation  may  be  accomplished  and 
there  Is  no  difference  If  doing  It  for  two  dimensional  or  for  axlsymmetrlc 


APPENDIX  A  -  THE  AXSYM  PROGRAM 


A.l  DIFFERENT  REGIONS  IN  THE  JET 

For  Che  two  dlaenelonel  Jet  with  Inlclel  Mach  number  greaCer  chan  unlcy 
Che  differenc  regions  are  ahoim  In  Figure  (16). 


Figure  16*  The  chree  regions  In  an  underexpanded  Jec. 

For  planar  2-0  flow  region  1  Is  a  uniform  flow  core,  region  2  is  a  simple 
region  in  which  only  one  family  of  characCeriscics  define  che  flow,  and  region 
3  which  concalns  Che  InCersecCion  of  Che  Cwo  families  of  characceclsCics. 
Because  our  incenclon  is  Co  find  soluCions  for  highly  underexpanded  Jets  wlch 
very  low  ambienc  pressure,  che  calculacion  of  furcher  downscream  flow  Is  noc 
necessary. 

For  che  axlsyimnecric  ring  jec,  we  use  che  same  definlclon  for  che 
differenc  regions  homver,  in  chis  case  none  of  che  chree  regions  has  uniform 
flow  and  is  noc  a  simple  region. 
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As  shown  in  Equations  (9,10)  Che  PM  function  ( v)  and  the  flow  direction 
(9)  of  a  point  at  location  l,j  may  be  calculated  from  the  v  and  6  of  two 
upstream  points  and  Later,  from  the  PM  function  at  the  new 

point  we  may  derive  the  local  Mach  number,  Che  local  pressure,  temperature, 
velocity  and  other  thermodynamic  parameters  as  required. 

Definition  of  the  mesh  of  points  for  the  different  regions  is  shown  in 
Figure  (17). 

A. 2  PROGRAM  FLOWCHART 

A  simplified  flowchart  for  Che  MOC  program  is  shown  in  Figure  (18).  The 
program  is  designed  Co  solve  both  axlsymmetric  as  well  as  two  dimensional  flow 
for  kD  »  2  it  solves  two  dimensional  flow 

for  kD  >  3  it  solves  axlsymmetric  flow  (This  is  also  the  default 
condition) 

Initial  data  such  as  ilach  number  and  pressure  at  the  exit  surface,  ambient 
pressure  and  jet  gas  parameters  are  input  data. 

Output  data  contains  the  following  for  each  mesh  point: 

Mach  number,  coordinates  of  mesh  point  (R,X),  flow  direction  ^,TETA)  , 
pressure,  temperature,  local  velocity,  Knudsen  number  basea  on  the 
distance  between  two  points  along  a  streamline,  mean  free  path  and 
breakdown  parameter  as  defined  by  Bird.* 

For  each  of  the  three  regions,  we  start  with  precalculated  boundary 
conditions  enabling  the  calculation  of  the  ilach  angles,  coordinates  of  raesh 
points  and  distances  d^  and  dri  as  described  in  [2]. 

*For  the  exit  plane  instead  the  Bird's  breakdown  parameter,  we  calculate  Che 
ratio  between  time  per  three  collisions  and  time  of  motion.  Sometimes  this 
ratio  may  be  regarded  as  a  measure  of  the  breakdown  of  the  continuum  theory. 


The  number  of  characteristics  used  in  the  program  is  arbitrary  and 
depends  on  the  required  resolution  (it  may  affect  also  the  accuracy  of  the 
results  and  the  amount  of  computation).  We  start  with  20  characteristics 
along  the  exit  cross  section  and  with  50  characteristics  in  the  Prandtl-Meyer 
fan  thus  a  total  of  70  characteristics  of  each  family  are  calculated.  In 
region  1  there  are  20  left  running  and  20  right  running  characteristics.  In 

region  2  there  are  20  right  running  and  50  left  running  characteristics.  In 

region  3  there  are  50  left  tunning  and  50  right  running  characteristics. 

In  region  3  we  limit  the  calculation  where  the  two  characteristics 

defining  a  new  mesh  point  Intersect  at  an  angle  smaller  than  the 

computational  accuracy.  In  fact  this  occurs  far  downstream  where  continuum 
theory  becomes  Invalid. 

After  defining  the  mesh  geometry  (successively)  we  calculate  the 
Prandtl-Meyer  function  and  flow  direction,  using  equations  (9,10). 

The  local  Mach  number  is  an  implicit  function  of  the  Pradtl-tleyer  angle. 
It  is  calculated  by  iterations  with  an  initial  guess  of  local  Mach  number 
set  equal  to  a  precalculated  Mach  number  at  an  adjacent  point,  and  the  slope 
of  the  function  as  given  by  Equation  (5).  Figure  (19)  shows  the  iterative 
procedure  for  evaluating  the  Mach  number  at  each  mesh  point. 


'^eSOL.T- 


Figure  19.  Iterative  procedure  for  Mach  number  calculation. 

Once  the  local  Mach  number  has  been  found  all  other  flow 
parameters  may  be  defined  using  the  Equation  (1,2,3). 


Figure  20.  The  oesh  of  characcerlsclcs  in  Region  2  and  Region  3 
AxisyMCric  ring  Jet.  Mg  -  4.  Altltude-200kni. 
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A.3  PROGRAM  'AXSYII'  LISTING 


AXSOOOlO 
AXS00020 
AXS00030 

THIS  PROGRAM  CALCULATES  THE  ISENTROPIC  EXPANSION  OF  A  JET  BY  MEANS  OFAXSOOOAO 
THE  METHOD  OF  CHARACTERISTICS. 

FOR  A  TMO  DIMENSIONAL  JET  'KD*  SHOULD  BE  SET  EQUAL  TO  2 

FOR  AN  AXISYMMETRIC  RING  JET  *KO*  SHOULD  BE  SET  EQUAL  TO  3 


$JOB 

C  PROGRAM  AXSYM 
C 
C 
C 
1C 
C 
iC 


AXS00050 
AXS00060 
AXS00070 
AXSOOOSO 

IMPLICIT  REAL»8(A-H,0-2,«)  AXS00090 
DIMENSION  TETA(20,SO),AN(20,50),R(20,50),XC20,50},PM(20,50>  AXSOOlOO 
DIMENSION  AMCOR(20.20),TETAC(20,20>,XC(20,20),RC(20,20),PMCC20,20)AXS00110 


DIMENSION  DENSF(20,50) 

DIMENSION  AMX(50,50),TETAX(50,50),XXC50,50),RX(50,50),PMX(50,50) 


AXS00120 

AXS00130 

AXS00140 

AXS00150 

AXS00160 

AXS00170 

AXS00180 

AXS00190 


TETA  IS  THE  FLOW  ANGLE  (RADIANS)  MEASURED  FROM  X  AXIS. 

AM  IS  THE  MACH  NUMBER 

R  IS  THE  RADIUS  (NORMAL  TO  THE  HALL) 

X  IS  THE  AXIAL  LOCATION  (PARALLEL  TO  THE  WALL) 

PM  IS  THE  PRANDTL  MEYER  FUNCTION 
CXXKXXXXXXXXXXKX1(XXXX)(XXXX)(XXXX)()(XKXKXXXXXK)(XX»XXXXXXXXXXXXKXXXKXXXXXXXXAXS00200 
C  THE  FOLLOWING  IS  DATA  FOR  THE  SPECIFIC  PROBLEM  AXS00210 

PAMB  =  8.4736E-5  AXS00220 

KD  =  3  AXS00230 

C  FOR  TWO  DIMENSIONAL  FLOW  KD=2  .-FOR  AXI SYMMETRICAL  FLOW  KD=3  AXS00240 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00250 


c 

AXS00260 

c 

CONSTANTS 

AXS00270 

PI  =  3.141593 

AXS00280 

BOLTZ  =  1.38032E-23 

AXS00290 

AVOG  »  6.0225E-i-26 

AXS00300 

RG  :  8314.3 

AXS00310 

c 

EXIT  SURFACE 

AXS00320 

AMO  =4.00 

AXS00330 

TO  =  300.0 

AXS00340 

PO  =  136.0 

AXS00350 

c 

AXS00360 

R1  =2.5 

AXS00370 

XI  =0.5 

AXS00380 

c 

AXS00390 

c 

GAS  DATA 

AXS00400 

GAMA  =  1.535 

AXS00410 

DIAM  =  2.95E-10 

AXS00420 

GM  =17.0 

AXS00430 

RJ  =  RG/GM 

AXS00440 

CXS  =  PIXDIAMXDIAM 

AXS00450 

GMM  =  GM/AVOG 

AXS00460 

c 

AXS00470 

c 

MESH  DEFINITION 

AXS00480 

c 

M1=CHARACTERISTICS 

FROM 

THE 

EXIT  PLANE 

AXS00490 

c 

N2=CHARACTERISTICS 

FROM 

THE 

CORNERS 

AXS00500 

N1  =  20 

AXS00510 

N2  =  50 

AXS00520 

c 

AXS00530 

c 

CONSTANTS  FOR  THE  ISENTROPIC  EXPANSION 

AXS00540 

A1  =  ( GAMA-1. 0)/GAMA 

AXS00550 

B1  =  1.0/Al 

AXS00560 

c 

AXS00S70 

A2  =DSQRT((GAMA+1.)/(GAMA-1.)) 

AXS00580 

B2  =  1.0/ A2 

AXS00590 

c 

AXS00600 

A3  =  (GAMA-1. )/2. 

AXS00610 

c 

AXS00620 

c 

C  =  MACHXMACHXA3  +  1. 

AXS00630 

c 

D  =  MACHXMACH  -1. 

AXS00640 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00650 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00660 


DEFINE  STAGNATION  PARAMETERS 


AXS00670 

AXS00680 

AXS00690 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  AXS00700 


C  s  (1.0-^A3XAM0XAM0) 
C  PRESSURE 


AXS00710 

AXS00720 


PN  >  CXXBl  AXS00730 

PST6  >  POXPN  AXS00740 

C  TEMPERATURE  AXS00750 

TSTO  «  TOXC  AXS00760 

C  DENSITY  AXS00770 

DO  *  PO/CRJXTO)  AXS00780 

'  DSTO  «  PSTG/(RJXTSTG)  AXS00790 

C  AXS00800 

CXXXJfXXXXXXXXX  AXS00810 

HRITE(6.1)PSTG,TST6,DST0.DIAM,GN  AXS00820 

1  FORMATC 'O', 'STAGNATION  PRESSURE** , E12. 5, '  TEMPERATURE*' ,F10. 5, AXS00830 

1'  DENSITY* ',E12. 5, '  MOL .DIAM*' , E12. 5, '  MOL .MASS*' ,F10.5)  AXS00840 

HRITE  (6,2)P0,T0,D0.AM0  AXS00850 

2  FORMATCO', 'EXIT  PLANE  PRESSURE* *,E12. 5, '  TEMPERATURE*', FIO. 5, AXS00860 

1'  DENSITY* ',E12. 5, '  MACH* * , FIO . 5)  AXS00870 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXAXS00880 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00890 
C  AXS0O9OO 

C  DEFINE  FREE  STREAM  PARAMETERS  AT  THE  RIGHT  CORNER  AXS00910 

C  AXS00920 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS00930 
C  MACH  NUMBER  AXS00940 

FSM  *  DSQRT(((PSTG/PAMB)XXA1-1.)/AS)  AXS00950 

C  TEMPERATURE  AXS00960 

FST  *  TSTG/(1.+A3XFSMXFSM)  AXS00970 

C  DENSITY  AXS00980 

FSD  *  PAMB/(RJXFST)  AXS00990 

C  MACH  ANGLES  FOR  HEAD  AND  TAIL  OF  FAN  AXSOIOOO 

AMIT  *DARSIN(l./AMO)  AXSOlOlO 

AMIH  =DARSIN(1./FSM)  AXS01020 

C  PRANDTL  MEYER  FUNCTION  FOR  HEAD  AND  TAIL  OF  FAN  AXS01030 

01  =DSQRT(AM0XAM0-1.)  AXS0I040 

D2  *DSQRT(FSMXFSM-1.)  AXS01050 

PMH  *  A2X0ATAN(B2XD2)-DATAN(D2)  AXS01060 

PMT  *  A2XDATAN(B2X01)-DATAN(D1)  AXS01070 

C  EXTERNAL  TURNING  ANGLE  (FREE  STREAM  ANGLE)=EXTA  AXS01080 

EXTA  *  PMH  -  PMT  AXS01090 

C  PRANDTL  MEYER  FAN  ANGLE  PMFA  AXSOllOO 

PMFA  *  EXTA  -  AMIH  ♦  AMIT  AXSOlllO 

C  AXS01120 

Cxxxxxxxxxxxxxx  AXS01130 

C  WRITE(6.3)PAMB,FST>FS0,FSM  AXS01140 

C  3  FORMATCO', 'FREE  STREAM  PRESSURE*', E12. 5, '  TEMPERATURE* ', FIO . 5, AXS01150 
C  1*  DENSITY*  *,E12. 5,  •  MACH*' ,  FI  0.5///)  AXSOlUO 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS0117O 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXSOllSO 
C  AXS01190 

C  DEFINE  THE  MESH  OF  CHARACTERISTICS-AND  FLOW  PARAMETERS  AXS01200 

C  AXS01210 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS01220 


1  FORMATCO', 'STAGNATION 
1'  DENSITY* ',E12. 5, ' 

WRITE  (6,2)P0,T0,D0,AM0 

2  FORMATCO', 'EXIT  PLANE 

1'  DENSITY* ',E12. 5,' 


PRESSURE**, E12. 5,' 
M0L.DIAM*',E12.5,' 


C  THE  CORNER  POINT  -  LEFT  RUNNING  CHAR. 

C  PMFA  IS  DEVIDED  INTO  N2  (*50)  NONEQUAL  DIVISIONS 
C 

RAT  *  (FSM/ AMO-1 )/0. 02 
RAT2  *  FL0AT(N2-1) 

El  *  0L0G(RAT)/DL0G(RAT2) 

C 

WRITE(6,4) 

4  FORMATCl', 'PRANDTL  MEYER  FAN  LINE  M/ 

1  P.M.  ANGLE  TETA'/) 

DO  20  N*2,N2 
EN2  *  FL0AT(N-2) 

ENl  *  FLOAT(N-l) 

AMB  *  AM0X(1.+.02X(EN2)XXE1) 

AMF  *  AM0X(1,+,02X(EN1)XXE1) 

DELM  *  AMF-AMB 
AM(1,N)  *  AMF 

D1  =DS0RT(AM(1,N)XAM(1,N)-1.) 

PM(1,N)  *  A2XDATAN(B2XD1)-DATAN(D1) 

AMI  *DARSIN(1./AM(1,N)) 

TETA(1,N)  *  PI/2.-(PM(l,N)-PMT) 

C  TETA  IS  THE  FLOW  ANGLE  ON  THE  CHARACTERISTICS  AT  THE  CORNER 


AXS01230 

AXS01240 

AXS01250 

AXS01260 

AXS01270 

AXS01280 

AXS01290 

AXS01300 

AXS01310 

AXS01320 

AXS01330 

AXS01340 

AXS01350 

AXS01360 

AXS01370 

AXS01380 

AXS01390 

AXS01400 

AXS01410 

AXS01420 

AXS01430 

AXS01440 


•*. 


ijyi  1.'“*  L"ii  ir*  l-*  v^ 


R(1,N)  *  R1  AXS01450 

X(1»N)  s  XI  AXS01460 

ALFAL  »  TETA(1,N)-^AMI  AXS01470 

C  ALFAL  IS  THE  ANCLE  OF  THE  LEFT  RUNNING  CHARACTERISTICS  AXS01480 

C  AXS01490 

C  PRESSURE.  TEMPERATURE  AND  DENSITY  VARIATION  AT  THE  CORNER  AXS01500 

C  =  AN(1,N)XAM(1,N)XA3-^1.  AXS01510 

PRES  >  PSTG/(CXXB1)  AXS01520 

TEMP  =  TSTG/C  AXS01530 

DENSF(l.N)  -  PRES/(RJ»TEMP)  AXS01540 

C  AXS01550 

CXXXXXXXX  AXS01560 

HRITE  (6.5)  N.AM(1.N).PM(1.N).TETA(1.N)  AXS01570 

5  FORMAT! •  •.1I5.3E20.5)  AXS01580 

CXXXXXXXK  AXS01590 

20  CONTINUE  AXS01600 

C  AXS01610 

CXXXXXXXXXXKXXXXXXXXXXXXXXXKXXKKXXXXKXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXAXS01620 
C  AXSD1630 

CXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS01640 
NRITE(6,10)  AXS01650 

10  FORMAT! >0'.'  I  J  MACH  R  X  TETA  AXS01660 

ITEMP  PRESS  VELOCITY  KNUDSEN  MFP  ND  AXS01670 

1  P')  AXS01680 

C  AXS01690 

C  AXS01700 

CxxxxxxxxxxxxxxxxxxxxxkxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS01710 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS01720 
C  AXS01730 

C  DEFINE  THE  PARAMETERS  AT  THE  EXIT  SURFACE  (PLANE)  AXS01740 

C  AXS017S0 

CXXXKXXXKXXXXXXXXXXXKXKKXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS  01760 


DO  25  J  =  l.Nl 
AMCOR! l.J)  =  AMO 
D  =  DSQRT!AM0XAM0-1. ) 

PMCd.J)  =  A2XDATAN!B2XD)-DATAN!D) 

25  TETACd.J)  =  PI/2. 

C 

C  THE  EXIT  PLANE  IS  DIVIDED  INTO  (Nl-1)  DIVISIONS  !N1  POINTS) 
DO  30  J3i,20 

SOUND  s  0SQRT!GAMAXRJKT0) 

VELO  3  AMC0R!1,J)XS0UND 
DNO  »  PO/!BOLTZXTO) 

FPO  =  .707/(DN0xCXS) 

DX  =  Xlx2./FL0AT(N1) 

AKNO  =  FPO/DX 
C 

XCd.J)  =  XlXd. -FLOAT!  J-1)/FL0AT!N1-1)X2.) 

CENTR  »  DABS(XC!1.J)) 

IF(CENTR.LT. 0.001)  XCd.J)  -  0. 

RCd.J)  =  R1 


C  KZ^lxxxxxxxxxxxxxxxxxxx 
IF  (J.GT.l)  GO  TO  29 

C  AT  THE  EXIT  PLANE  THE  BREAKDOWN  PARAMETER  IS  EVALUATED  BY  MEANS  OF 
C  THE  RATIO  OF  THE  COLLISION  TIME  AND  THE  FLOW  TIME.  LENGTH  SCALE  IS 
C  THE  MESH  DIMENSION. 

SCALE  3  OXXDTAN(DARSIN!l./AMO))/2. 

TIMEl  =  SCALE/VELO 

TIME2  s  3./(4.XCXSXDNOXDSORT!BOLT2XTO/!PlxGMM))) 

P  =TIME2/TIME1 
DENSFCl.l)  -  DO 
29  CONTINUE 


AXS01770 

AXS01780 

AXS01790 

AXS01800 

AXS01810 

AXS01820 

AXS01830 

AXS01840 

AXS01850 

AXS01860 

AXS01S70 

AXS01880 

AXS01890 

AXS01900 

AXS01910 

AXS01920 

AXS01930 

AXS01940 

AXS01950 

AXS01960 

AXS01970 

AXS01980 

AXS01990 

AXS02000 

AXS02010 

AXS0202Q 

AXS02030 

AXS02040 

AXS020SO 

AXS02060 

AXS02070 

AXS02080 


C  AXS02090 

WRITE!6. 11 )1.J.AMC0R!1.J).R1.XC! l.J). TETACd.J), TO. PO. VELO, AKNO. FPAXS02100 
10, DNO, P  AXS02110 

11  FORMAT!'  •.2I4,5F10.3,3E12.3,F9.4,2E12.3)  AXS02120 

30  CONTINUE  AXS02130 

C  AXS02140 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXKXXXXXXXXXXXAXS02150 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXKXXXXXXKXXAXS02160 


I 


non  oooo  o  ooo  o  o  o  o  oooooo 


AXS02170 
AXS02180 
AXS02190 
AXS02200 
AXS02210 

KXXXX»XX)fXlli(KXXXXXXi(XXXXXXXXXXXX]fXXXXillflfXXXXXXXXXXXl(XllXlCXXXifXXXXXXXKXXXAXS02220 


CALCULATE  THE  FLOW  PARAMETERS  IN  THE  CORE  BOUNDED  BY  THE  THO  MACH 
WAVES  STARTING  AT  THE  NOZZLE  LIPS  (CORNER  POINTS) 

AT  THE  EXIT  THE  FLOW  IS  ASSUMED  TO  BE  UNIFORM 


198 


WRITE  (6,198) 
FORMATCl',  • 
APMsO. 

BPM-0. 

CPM>0. 


CORE 


•/) 


AXS02230 
AXS02240 
AXS0225a 
AXS02260 
AXS02270 
AXS02280 
AXS02Z90 
AXS02300 
AXS02310 
AXS02320 
AXS02330 
AXS02340 
AXS02350 
AXS02360 
AXS02370 
AXS02380 

XC(I,J)  =  (RC(I-1, J)-RC(I-1,J+1)+XC(I-1,J)XDTAN(ALFAL)+XC(1-1, J+1)AXS02390 


DO  199  I  >  2,N1 
WRITE  (6,10) 

DO  199  J  >  1,N1 
IF  ((I+J-1).GT.N1)  GO  TO  199 
AMIL  >  DARSIN(1./AMC0R(I-1,J)) 
ALFAL  =  PI-(TETAC(I-1,J)+AMIL) 

AMIR  -  0ARSIN(1./AMC0R(I-1,J4'1)) 
ALFAR  =  TETAC(I-1,J+1)-AMIR 


1XDTAN( ALFAR))/(DTAN( ALFAL )+DTAN(ALFAR)) 

CENTR  =  DABS(XC(I, J)) 

IF  (CENTR.lt. 0.001)  XC(I,J)  *  0. 

RC(I,J)  =  RC(1-1,J+1)+(XC(I,J)-XC(I-1, J+1))XDTAN(ALFAR) 

DKSI  =  DSQRT((XC(I,J)-XC(I-l,J+l))«3f2+(RC(I,J)-RC(I-l,J+l))xx2) 
DETA  =  DS()RT((XC(I,J)-XC(I-1,J))«2+(RC(I,J)-RC(I-1,  J))XX2) 

CALCULATE  NOW  THE  PRANOTL  MEYER  FUNCTION  AND  FLOW  ANGLE  IN  CORE. 

APM  =  PMC(I-1, J)+PMC(I-1,J+1)+TETAC(I-1, J+1)-TETAC(I-1,J) 

IF  (KD.EQ.2)  GO  TO  151 

BPM  =  DSIN(AMIR)XDSIN(TETAC(I-1,J+1))/RC(I-1,J+1)»DKSI 
CPM  =  DSIN(AMIL)»DSIN(TETAC(I-1,J))/RC(I-1,J)XDETA 
PMC(1,J)  =  (APM+BPM+CPM)/2.0 

APM  =  PMC(I-1,J+1)-PMC(I-1,J)+TETAC(I-1,J+1)+TETAC(I-1,J) 
TETAC(I,J)  =  (APM+BPM~CPM)/2.0 


151 


DEFINE  MACH  NUMBER  (BY  ITERATIONS) 
INITIAL  GUESS  AMCORd, J)=AMCOR(I-l, J+1) 

AMG  =  AMCOR(I-l,J+l) 

KZ  =  0 

154  IF  (KZ.GE.lOO)  GO  TO  160 
KZ  =  KZ+1 
C  *  AMGXAMGXA3+1 . 

D  =  DSQRT(AMGXAMG-1.) 

PMCAL  *  A2XDATAN(B2XD)-DATAN(D) 
DELNI  =  PMCAL  -  PMC(I,J) 

DEL  -  DABS( DELNI) 

IF  (DEL. LT. .000002)  GO  TO  160 
IF  (DELNI. LT. 0. )  GO  TO  156 
AMG  =  AMGX.999 
GO  TO  154 

156  AMG  =  AMGX(1.-DELNIXC/D) 

GO  TO  154 

160  AMCOR(I,J)  =  AMG 


CALCULATE  FLOW  PARAMETERS 
IF  (J.GT.l)  GO  TO  197 
C  =  AMC0R(I,J)XAMC0R(I,J)XA3-M. 
PRES  =  PSTG/(CxxBl) 

TEMP  =  TSTG/C 

DN  =  PRES/(BOLTZXTEMP) 

FP  =  .707/(DNXCXS) 

SCALE  =  DKSI  X  DSIN(ALFAL) 

AKN  =  FP/SCALE 


AXS02400 

AXS02410 

AXS02420 

AXS02430 

AXS02440 

AXS02450 

AXS02460 

AXS02470 

AXS02480 

AXS02490 

AXS02500 

AXS02510 

AXS02520 

AXS02530 

AXS02540 

AXS02550 

AXS02560 

AXS02570 

AXS02580 

AXS02590 

AXS02600 

AXS02610 

AXS02620 

AXS02630 

AXS02640 

AXS02650 

AXS02660 

AXS02670 

AXS02680 

AXS02690 

AXS02700 

AXS02710 

AXS02720 

AXS02730 

AXS02740 

AXS02750 

AXS02760 

AX302770 

.AXS02780 

AXS02790 

AXS02800 

AXS02810 

AXS02870 

AXS02830 

AXS02840 

AXS02850 

AXS02860 

AXS02870 

AXS02880 


SOUND  -  DSQRT(6AMAXRJXTEMP)  AXS02890 

VEL  s  SOUNOxANCOR(I,J)  AXS02900 

C  AXS02910 

C  BREAKDOWN  PARAMETER  AS  DEFINED  BY  'BIRD*.  AXS02920 

DENSF(I,J)  *  PRES/(RJ»TEMP)  AXS02930 

DDENS  =  DENSF(I-1,J)  -  DENSF(I.J)  AXS02940 

COLF  s  4.XCXSXDNXDSQRT(B0LTZXTENP/(PIXGMM))  AXS02950 

P  =  VELXDDENS/(SCALEXDENSF(I,J)XCOLF)  AXS02960 

197  CONTINUE  AXS02970 

C  AXS02980 

C  TIMEl  =  SCALE/VEL  AXS02990 

C  TIME2  -  3./(4.XCXSXDNXDSQRT(B0LTZXTEMP/(PIXGMM)))  AXS03000 

C  P  =  TIME2/TIME1  AXS03010 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS03020 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03030 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS 03040 
WRITE(6, 11)1. J, AMCORCI, J), RCCl, J), XCd, J),TETAC(I,J), TEMP, PRES, VELAXS03050 
1,AKN,FP,DN,P  AXS03060 

199  CONTINUE  AXSi)3070 

C  AXS03080 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS03090 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS03100 
C  MATCH  CORE  AND  FAN  POINTS  AXS03110 

DO  200  1=1, N1  AXS03120 

X(I,1)  =XC(I,1)  AXS03130 

R(I,1)=  RC(I,1)  AXS03140 

AM(I,1)  =  AMCOR(I,l)  AXS03150 

PM(I,1)  =  PMC(I,1)  AXS03160 

200  TETA(I,1)  =  TETACd,!)  AXS0S170 

C  AXS03180 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAXS03190 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03200 
C  AXS03210 

C  CALCULATE  FLOW  PARAMETERS  IN  REGION  2  (SIMPLE  PRANDTL  MEYER  FAN).  AXS03220 

C  AXS03230 

CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS03240 


WRITE  (6,298)  AXS03250 

298  FORMAT  Cl','  REGION  2*/)  AXS03260 

AXS03270 

KFIN  =  51  AXS03280 

AXS03290 

DO  299  I  =  2,N1  AXS03300 

WRITE(6,10)  AXS03310 

DO  299  J  =  2,N2  AXS03320 

IF  (J.GT.KFIN)  GO  TO  299  AXS03330 

KZ  =  0  AXS03340 

AMIL  =DARSIN(1./AM(I-1, J))  AXS03350 

AMIR  =DARSIN(1./AM(I,J-1))  AXS03360 

ALFAL  =  PI-(TETA(I-1, J)+AMIL)  AXS03370 

ALFAR  =  TETAd,  J-1)-AMIR  AXS03380 

AXS03390 

CHECK  ANGLES  AND  CALCULATE  COORDINATES  AXS03400 

ANGLEl  =  PI/2. -.000001  AXS03410 

ANGLE2  =  PI/2. -t-. 000001  AXS03420 

IF  (ALFAL. LE. ANGLEl. OR. ALFAL. GE.ANGLE2)  GO  TO  201  AXS03430 

X(I,J)  =  X(I-1,J)  AXS03440 

GO  TO  207  AXS03450 

201  IF  (ALFAR. LE. ANGLEl. OR. ALFAR. GE.ANGLE2)  GO  TO  205  AXS03460 

X(I,J)  =  X(I,J-1)  AXS03470 

GO  TO  207  AXS03480 

AXS03490 

205  Xd,J)  =  (Rd-l,  J)-R(I,J-l)+X(I,J-l)XDTAN(ALFAR)+Xd-l,  J)XDTAN(ALFAXS03500 
lAL ))/( DTAN( ALFAL )+DTAN( ALFAR))  AXS03510 


REGION 


DO  299  I  =  2,N1 
WRITE(6,10) 

DO  299  J  =  2,N2 
IF  (J.GT.KFIN)  GO  TO  299 
KZ  =  0 

AMIL  =DARSIN(1./AM(I-1, J)) 

AMIR  =DARSIN(1./AM(I,J-1)) 

ALFAL  =  PI-(TETA(I-1, J)+AMIL) 

ALFAR  =  TETAd,  J-1)-AMIR 
C 

C  CHECK  ANGLES  AND  CALCULATE  COORDINATES 
ANGLEl  =  PI/2. -.000001 
ANGLE2  =  PI/2. -t-. 000001 

IF  (ALFAL. LE. ANGLEl. OR. ALFAL. GE.ANGLE2)  GO  TO  201 
X(I,J)  =  X(I-1,J) 

GO  TO  207 

201  IF  (ALFAR. LE. ANGLEl. OR. ALFAR. GE.ANGLE2)  GO  TO  205 
X(I,J)  =  X(I,J-1) 

GO  TO  207 


207  IF  (ALFAL. LE. 0.000001. OR. ALFAL. GE. 0.000001)  GO  TO  209 
R(I,J)  =  R(I-1,J) 

GO  TO  213 

209  IF  (ALFAR. LE. 0.000001. OR. ALFAR. GE. 0.000001)  GO  TO  211 
R(I,J)  =  R(I,J-1) 

GO  TO  213 

211  R(I,J)  =  (X(I.J)-X(I, J-1))XDTAN(ALFAR)+R(I,J-1) 


AXS03520 

AXS03530 

AXS03540 

AXS03550 

AXS03560 

AXS03570 

AXS035S0 

AXS03590 

AXS03600 


oo  oo  o  oo  on  o  o  ooo  o  oooo 


213  DKSI  »  DSQRT((R(I,J)-R(I,J-1))»*2+(XCI,J)-X(I,J-1))K«2) 

DETA  »  DSQRT((R(I,J)-R(I-1,J))*K2+<XCI,J)-X(I-1,J))*K2) 

IF  (R(I,J).6T.O..AND.X(I,J).GT.O.)  GO  TO  219 

KFIN  =  J-1 

HRITE(6,12) 

12  FORMAT! •  FURTHER  POINTS  ON  THE  CHARACTERISTICS  ARE  DIVERGENT* 
219  CONTINUE 

LOCATION  OF  THE  NEH  MESH  POINT  HAS  BEEN  FOUND 

CALCULATE  NON  PRANDTL  MEYER  FUNCTION  AND  FLOW  DIRECTION  FOR  NEW  POINT 
APM  =  PMCI,J-1)+PMCI-1,J)-TETACI-1,J)+TETA(I,J-1) 

IF  (KD.Eq.2)  GO  TO  251 

BPM  »  DSIM(AMIR)*DSIN(TETA(I,J-1))/R(I, J-l)3fDKSI 
CPM  =  DSIN(AMIL)»DSIN(TETA(I-1,J))/R(I-1,J)XDETA 
251  PM(I,J)  =  .5XCAPM+BPM+CPM) 

APM  =  PM(I, J-1)-PM(I-1, J)+TETACI,J-1)+TETA(I-1, J) 

TETA(I,J)  =  .5*(APM+BPM-CPM) 

CALCULATE  NOW  MACH  NUMBER  FOR  EACH  POINT 
INITIAL  GUESS  AMd.J)  =  AM(I-1,J) 

AMG  =  AM(I-1,J) 

KZ  =  0 

254  IF  (AMG.GT.200.)  GO  TO  257 
D  =  DS0RT(AMGXAMG-1. ) 

GO  TO  258 

257  D  =  AMG 

258  IF  (KZ.GE.lOO)  GO  TO  260 
KZ  =  KZ  +  1 

PMCAL  =  A2XDATAN<B2KD)-DATAN(D) 

DELNI  =  PMCAL  -  PMCI.J) 

DEL  =  DABS! DELNI) 

IF  !DEL.LT. .000002)  GO  TO  260 
IF  !DELNI.LT.O. )GO  TO  256 
AMG  =  AMG*. 999 
GO  TO  254 
'/////////* 

256  IF  !AMG.LT.2000.)GO  TO  2560 
KFIN  *  J-1 
GO  TO  299 


AXS03610 

AXS03620 

AXS03630 

AXS03640 

AXS03650 

AXS03660 

)AXS03670 

AXS03680 

AXS03690 

AXS03700 

AXS03710 

AXS03720 

AXS03730 

AXS03740 

AXS03750 

AXS03760 

AXS03770 

AXSJ)3780 

AXS03790 

AXS03800 

AXS03810 

AXS03820 

AXS03830 

AXS03840 

AXS03850 

AXS03860 

AXS03870 

AXS03880 

AXS03890 

AXS03900 

AXS03910 

AXS03920 

AXS03930 

AXS03940 

AXS03950 

AXS03960 

AXS03970 

AXS03980 

AXS03990 

AXS04000 

AXS04010 

AXS04020 


/////////////////////////////////////////////////////////////////////// AXS04030 


2560  D1  =  !A3*AMGXAMG+1. ) 

AMG  =  AMGX!1.-DELNI*D1/D) 

GO  TO  254 

260  AM!I,J)  =  AMG 

CALCULATE  NOW  LOCAL  TEMPERATURE, PRESSURE, VELOCITY, KNUDSEN  NO. 

C  =  AM!I,J)XAM!I, J)XA3+1 
PRES  =  PSTG/!CxxB1) 

TEMP  =  TSTG/C 

DN  =  PRES/!BOLTZXTEMP) 

DENSF!I,J)  =  PRES/!RJXTEMP) 

DDENS  =  DENSF!I-l, J-1)-DENSF!I, J) 

SCALE  =  DS0RT!!X!I,J)-X!I-1,J-1))xx2+!R!I,J)-R!I-1,J-1))**2) 
XXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXKKXXXXXXXXXX 

271  FP  =  .707/!DNXCXS) 

AKN  =  FP/SCALE 

SOUND  =DSQRT!GAMAXRJXTEMP) 

VEL  =  AM!I,J)XS0UND 


COLF  =  4.xCXSXDNXDS0RT!B0LTZxTEMP/!PIXGMM)) 

P  =  VELXDDENS/!SCALExDENSF!I, J)xC0LF) 

PRINT  RESULTS  FOR  MESH  POINTS 
KL  =  !-l)xxl 
IF  !KL.LT.0)  GO  TO  299 

WRITE!6,11)I, J,AM!I, J),R!I,J),X!I,J),TETA!I,J),TEMP,PRES,VEL,AKN, 


AXS04040 

AXS04050 

AXS04060 

AXS04070 

AXS04080 

AXS04090 

AXS04100 

AXS04110 

AXS04120 

AXS04130 

AXS04140 

AXS04150 

AXS04160 

AXS04170 

AXS04180 

AXS04190 

AXS04200 

AXS04210 

AXS04220 

AXS04230 

AXS04240 

AXS04250 

AXS04260 

AXS04270 

AXS04280 

AXS04290 

AXS04300 

AXS04310 

AXS04320 


1FP,DN,P  AXS04330 

299  CONTINUE  AXS04340 
C  AXS04350 
CXXXi()iXXXi(X)fXXXX]fXXXXXXXKXXXIfXX]fXXXXI(Xl(XXXXI(KXKX](XXXXXXXXXXXXXKXXXXXXKKXAXS04360 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxkxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxaXS04370 
C  AXS04380 
C  NATCH  'REGION  2*  AND  'REGION  3'  POINTS  AXS04390 
C  AXS04400 
CXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXKXXXXKXKXXXXXXKXXKXXKXXXXXXXKXXXXXXXXAXS04410 

L  -  KFIN-  1  AXS04420 
DO  300  J  s  1,L  AXS04430 
XX(1,J)  s  X(20,J)  AXS04440 
RX(1,J)  =  R(20,J)  AXS04450 
AMX(1,J)  -  AM(20,J)  AXS04460 
PMX(1,J)  -  PM(20.J)  AXS04470 

300  TETAXCl.J)  »  TETAC20,J)  AXS04480 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxKxxxxxxxxxxxAXS 04490 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS04500 

C  AXSD4510 
C  CACULATE  FLOH  PARAMETERS  FOR  REGION  3  AXS04520 
C  AXS04530 
CxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxAXS04540 


WRITE  (6.397) 

397  FORMAT  Cl','  REGION  3'/) 

DO  399  1  -  2,L 
DO  399  J  =  I,L 
IF  (J.GT.KFIN)  GO  TO  399 
IF  (J.GT.I)GO  TO  320 
WRITE  (6.10) 

320  KZ  =  0 

AMIL  =  DARSINd./AMXd-l.  J)) 

ALFAL  =  PI-(TETAX(I-1. J)+AMIL) 

IF  (J.GT.I)  GO  TO  301 
ALFAR  =  ALFAL 
TETAXd.J)  =  PIX.5 
TETAXd.J-l)  =  PI-TETAX(I-1,J) 

XXd.J)  =  0. 

RXd.J)  =  RX(I-1,J)+XX(I-1,  J)XDTAN(ALFAL) 
RX(I,J-1)  s  RX(I-1,J) 

XX(I,J-1)  =  -XXd-l.J) 

DKSI  =  (RX(I,J)-RX(I-1,J))/DSIN(ALFAL) 

DETA  =  DKSI 

PMXd.J-l)  =  PMX(I-1,J) 

GO  TO  316 

301  AMIR  =  DARSINd./AMXd,  J-D) 

ALFAR  =  TETAX(I.J-1)-AMIR 

IF  (ALFAL. LE.ANGLEl .OR. ALFAL. GE.ANGLE2)  GO  TO  302 
XXd.J)  *  XXd-l.J) 

GO  TO  307 

302  IF(ALFAR. LE.ANGLEl. OR. ALFAR. GE.ANGLE2)  GO  TO  305 
XXd.J)  *  XX(I,J-1) 

GO  TO  307 


AXS04550 

AXS04560 

AXS04570 

AXS04580 

AXS04590 

AXS04600 

AXS04610 

AXS04620 

AXS04630 

AXS04640 

AXS04650 

AXS04660 

AXS04670 

AXS04680 

AXS04690 

AXS04700 

AXS04710 

AXS04720 

AXS04730 

AXS04740 

AXS04750 

AXS04760 

AXS04770 

AXS04780 

AXS04790 

AXS04800 

AXS04810 

AXS04820 

AXS04830 

AXS04840 

AXS04850 

AXS04860 


305  XXd.J)  =  (RXd-l. J)-RX(I, J-1)+XX(I,J-1)XDTAN(ALFAR)+XX(I-1,J)XDTAAXS04870 


1N(  ALFAL  ))/(DTAN(  ALFAL  )-i-DTAN(  ALFAR)) 

307  IF  (ALFAL. LE. 0.000001. OR. ALFAL. GE. 0.000001)  GO  TO  309 
RXd.J)  =  RX(I-1,J) 

GO  TO  315 

309  IF  (ALFAR. LE. 0.000001. OR. ALFAR. GE. 0.000001)  GO  TO  311 
RXd.J)  =  RXd.J-l) 

GO  TO  315 

311  RXd.J)  =  RX(I,J-1)  +  (XX(I,J)-XX(I,J-1))XDTAN(ALFAR) 

315  DKSI  =  DS0RT((RX(I, J)-RX(I,J-1))x»2*(XX(I. J)-XX(I. J-1))XX2) 
DETA  =  DSQRT((RX(I,J)-RX(I-1,J))Xx2+(XX(I.J)-XX(I-1.J))Xk2) 

316  CONTINUE 

IF  (RX(I.J).GE.O.  .AND. XXd.J). GE.O.)GO  TO  319 
KFIN  =  J-1 
WRITE  (6.398) 


AXS04880 

AXS04890 

AXS04900 

AXS04910 

AXS04920 

AXS04930 

AXS04940 

AXS04950 

AXS04960 

AXS04970 

AXS04980 

AXS04990 

AXS05000 

AXS05010 

AXS05020 

AXS05030 

AXS05040 


398  FORMAT  (•  S 'FURTHER  POINTS  ON  THE  CHARACTERISTICS  ARE  DIVERGENT* 
C 

C  LOCATION  OF  THE  NEW  POINT  HAS  BEEN  FOUND 

CXXXXKXXXXXXKXXXXKXKXXXXXXXXXXKXXXKXXXKXXlfXXXXXKXXXXXXKXXXXXXXXKXXXXXX 

C 

C  CALCULATE  NOW  P.M.  ANGLE  AND  FLOW  DIRECTION 

319  APM  =  PMX(I,J-1)+PMX(I-1,J)-TETAX(I-1,J)+TETAX(I,J-1) 

IF(KD.EQ.2)  GO  TO  351 

BPM  =  DSIN(AMIR)xDSIN(TETAX(I,J-l))/RX(I,J-l)xDKSI 
CPM  =  DSIN(AMIL)XDSIN(TETAX(I-l,J))/RXtI-l, J)KDETA 
351  PMXd.J)  =  .5X(APM+BPM+CPM) 

APM  »  PMXCI,J-1)-PMX(I-1,J)+TETAX(I,J-1)+TETAX(I-1,J) 

TETAXd.J)  =  .5*(APM+BPM+CPM) 

C 

C  CALCULATE  NON  THE  MACH  NUMBER  FOR  EACH  POINT 
C  INITIAL  GUESS  AMXCI.J)  =  AMXd-l,J) 

AMG  -  AMXd-l»J) 

KZ  -  0 

KH=0 

KL=0 

354  IF  (AMG. GT. 2000.)  GO  TO  357 
D  =  DSQRT(AMGXAMG-1.) 

GO  TO  358 

357  D  =  AMG 

358  IF  (KZ.GE.50)GO  TO  360 
KZ  =  KZ  +1 

PMCAL  =  A2XDATAN(B2»D)-DATAN(D) 

DELNI  =  PMCAL-PMXd,J) 

DEL  =  DABS( DELNI) 

IF  (DEL. LT. .000002)  GO  TO  360 
IF  (DELNI. LT. 0. )GO  TO  356 
AMG  =  AMGX.98 
GO  TO  354 

356  IF  (AMG. LT. 5000.)  GO  TO  3560 
KFIN  =  J-1 
GO  TO  399 

3560  D1  =  A3XAMGXAMG+1 . 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DDELNI  =  DELNI 

IF  (AMG.GT.20.)  DDELNI  *  DELNIX( .95»KKZ) 
Cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
AMG  s  AMGXd.-DDELNiXDl/D) 

GO  TO  354 

360  AMX(I,J)  »  AMG 
C 

C  CALCULATE  NOW  THE  LOCAL  TEMPERATURE,  PRESSURE,  VELOCITY, KNUDSEN 
C  =  AMX(I,J)*AMX(I,J)XA3+1. 

PRES  =  PSTG/(CXXB1) 

TEMP  =  TSTG/C 

DN  =  PRES/(BOLTZ*TEMP) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DENSF(I,J)=PRES/(RJXTEMP) 

DDENS=DENSF(I-1, J-1)-DENSF(I,J) 

FP  =  .707/(DNXCXS) 

SCALE=DSORT((XX(I, J)-XX(I-1,J-1))**2+(RX(I, J)-RX(I-1, J-1))XX2) 
AKN=FP/SCALE 
C 
C 

CXXXXXXXXXXXXXXXXKXXXXXXKXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXKXXXXXXXX 

C  ALFALL  =  -ALFAL+PI 

C  AF  =  ALFALL-PIy'2. 

C  AF  =  DABS(AF) 

C  IF  (AF.LT. .000001)  GO  TO  370 

C  BF  s  l./DSORTd.  +  DTAN(ALFALL)»«2) 

C  SCALE=BF*(RX(I,J-1)-DTAN(ALFALL)X(XX(I, J-1)-XX(I-1, J))-RX(I-1, J)) 

.C  SCALE  -  DABS(SCALE) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXX 

SOUND  =  DSQRT(GAMAXRJXTEMP) 

VEL  =  AMX(I,J)XS0UND 
C 

C0LF=4 . XCXSXDNXDSORK  B0LTZ»TEMP/(PIXGMM) ) 
P=VELXDDENS/(SCALEXDENSF(I, J)XCOLF) 
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AXS05120 

AXSOSISO 

AXS05140 
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AXS05520 

AXS05530 

AXS05540 

AXS05550 

AXS05560 

AXS05570 

AXS05580 

AXS05590 

AXS05600 

AXS05610 

AXS05620 

AXS05630 

AXS05640 

AXS05650 

AXS05660 

AXS05670 

AXS05680 

AXS05690 

AXS05700 
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AXS05720 

AXS05730 

AXS05740 

AXS05750 

AXS05760 


I 

•Sy 


'  w"* 

*  V  • 
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> 


u  ,v' 


or 


PRINT  RESULTS 

MRITE<6,  11)1,J.AMX(I,J),RX(I,J),XX(I 
1AKN,FP,DM,P 
399  CONTINUE 
STOP 
END 

GENTRY 


AXS05770 

AXS05780 

, J),TETAX(I,J), TEMP, PRES, VEL,AXS05790 

AXS05800 

AXS05810 

AXS05820 

AXS05830 

AXS05840 


A. 4  LIST  OF  MAIN  SYMBOLS  USED  IN  'AXSYM* 


Physical  Name 


ambient  pressure 


K  dimensions 


Units 


pascals  I  real 


Description 


ambient  atmosphere 


Boltzmann  constant  I  Joule 


integer 


real 

constant 


real 


Avogadro's  constant  molecules 

mol 


Universal  gas 
constant 


Mach  No.  (E.P) 


Temperature  (E.P) 


Joule 


mol.de 


Pressure  (E.P) 


Cylinder  radius 


ressure 


6.0225  X  10 
1/kmol 


8314.3 


Mach  number  at  exit 
plane 


Temperature  at  exist 
plane 


Pressure  at  exit 
lane 


Ra 


Half  width  of  nozzle 


averaged  heat 
capacity  ratio  (let) 


averaged  molecular 
diameter  (jet) 


averaged  molecular  | 
weight  (jet) 


gas  constant  (jet) 


collision  cross 
section  (hard  sphere) 


averaged  mass  of  a 
molecule 


number  of  divisions 
(characteristics  from 
the  exist  plane) 


number  of  character¬ 
istics  in  the  Prandtl 
Meyer  fan 


A. 4  LIST  OF  IIAIN  SYMBOLS  USED  IN  'AXSYM'  (CONTINUED) 


Parameter 

Name 

Physical  Name 

Units 

Type 

Description 

A2 

real 

l( Y  +  l)/( Y  -  1)1 

B2 

real 

1/A2 

A3 

real 

Y-  1 

2 

C 

real 

M^  Y  -  1  +  1 

2 

D,D1  ,D2 

real 

m2  -  1 

PSTG 

real 

stagnation  pressure 

TSTG 

“k 

real 

stagnation  tempera¬ 
ture 

DSTG 

real 

stagnation  density 

FSM 

real 

free  stream  Mach 
number 

FST 

“k 

_ 

real 

free  stream  tempera¬ 
ture 

FSD 

real 

free  stream  density 

AMIT 

MT 

radius 

real 

Mach  angle  (tail  of 
P.M.  fan) 

AMIH 

UH 

radius 

real 

Mach  angle  (head  of 
P.M.  fan)  i 

PMT 

Vt 

radius 

real 

P.M.  function  (tail)  | 

PMH 

radius 

real 

P.M.  function  (head)  | 

KAT,RAT2, 

El 

real 

used  to  define  a  | 

loparithmic  division 
of  the  corner  char-  | 
acteristics  (50  lines] 
in  P.M.  fan)(a  linear] 
division  would  have  | 
resulted  in  concen-  | 
tration  of  character- | 
is  tics  at  high  tiach  | 
numbers)  ] 

UELM 

real 

difference  of  Mach  | 

numbers  | 

Ail(  I,J) 

M 

real  2-D 
array 

Mach  number  at 
location  (l,j)(used  | 
for  the  corner) 

PM(I,J) 

V 

radius 

P.M.  function  at 
location  (l,j)(used 
for  the  corner) 

ami 

M 

radius 

real 

Mach  angle 

R(I,J) 

m 

radius  (or  ordinate) 
at  (I,j) 

real  2-D | flow  direction  at 
array  I  (I,i) 


r—' 


TETA( I,J) 


radius 


A. 5  AXSYM  PROGRAM  USER'S  GUIDE 


1.  Input  data: 

ambient  pressure  PAMB 

Mach  number  (Exit  plane)  AMO 

Temperature  (Exit  plane)  TO 

Pressure  (Exit  plane)  PO 

Half  width  of  nozzle  K1 

Radius  of  nozzle  ring  R1 

Specific  heat  ratio  of  jet  gas  GAMA 

average  molecular  diameter  (jet)  DIAM 

average  molecular  weight  GM 

2.  Options  for  flow  geometry 

two  dimensional  flow  KD~2 

ring  jet  KD“3 

default  condition  KD*3 

3.  Resolution  of  mesh  points 


to  change  the  resolution  of  the  mesh  points 
a  -  change  N1  and  N2  as  necessary 

b  -  change  'DIMENSIONS'  according  to  new  values  of  Nj  and  N2 
c  -  define  distribution  of  ilach  lines  in  the  Prandtl-tleyer  fan  as  required 
(program  lines  126-131) 

4.  Execution  commands: 

After  copying  program  into  USER'S  FILE: 

WATFIV  AXSYM  *  (XTYPE 

The  program  will  run  on  user's  terminal  under  WATFIV.  A  soft  copy  of 
program  listing  and  output  listing  will  be  stored  in  user's  disk  named 
AXSYM  LISTING. 
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Hard  Copy 

PRINT  AXSYM  LISTING. 

Program  Outputs. 

All  necessary  outputs  are  automatically  listed  by  the  program. 

Figure  (20)  and  Figure  (21)  show  the  resulting  mesh  of  characteristics 
calculated  for  an  altitude  of  200  km  for  110=4  and  H0=2  respectively. 

In  these  figures  we  also  show  some  Isotherms  and  the  limit  where  the 
breakdown  parameter  equals  0.05.  These  lines  are  plotted  (manually)  using 
interpolation  procedures.  Data  along  the  breakdown  line  is  input  data  for 


the  molecular  flow 


APPENDIX  B  SIMUL  PROGKAtl 


B.l  DATA  ORGANIZATION 

Because  of  the  large  number  of  molecules,  cells,  regions  and  sectors  In 
the  simulation  and  the  large  number  of  data  related  to  each  molecule,  each 
cell  and  region  to  be  stored,  special  precautions  should  be  taken  in  order  not 
to  overflow  the  available  computer  memory. 

The  following  data  organization  was  used  in  SIHUL.  Figure  22  shows  the 
geometry  of  one  sector. 


Figure  22.  Definition  of  sector  geometry  and  cell  volume. 

cell  volume  -  v(I)  -  R(I)  *  DDALFA  x  DR  x  RSI  *  DFI 

v(I)  _  R(I)  »  RSI(I) 
v(l)  “  R(l)  •  RSI 


v(l}  is  the  volume  of  smallest  cell  in  a  region 


Tables  of  Molecules  and  Their  Parameters 


P1(L,N1M)  -  Light  molecules  of  the  jet 

P2(L,N2M)  -  Heavy  molecules  of  the  jet 

P3(L,N3M)  -  Ambient  gas  molecules  (not  used  in  the  present  program) 

NIM,  N2M,  N3M  -  maximum  number  of  molecules  in  simulation.  Number  of 
active  molecules  may  be  smaller  or  equal  to  (N<I>M). 
L=l,2,3  -  cartesian  components  of  velocity  v^,  Vy,  V2  [m/s] 

L=4  -  radial  coordinate  [meters] 

if  r=-99  the  particular  molecule  is  inactive 
L»5  -  angular  coordinate  [radians]. 

This  table  is  generated  each  time  the  simulation  is  initiated  for  a  region. 
That  means,  the  same  group  of  molecules  (as  stored  in  the  computer  memory)  is 


used  to  simulate  the  flow  in  all  regions  in  the  computation  domain. 


IP(N1A+N2A+N3A)  Table  of  the  addresses  of  the  active  molecules 
arranged  la  order  of  their  species  and  In  the  order  of  their  cells 
IC(N,I,j)  -  table  of  cells  (in  a  region)  integer  data 
I  =  1,10  -  radial  index 

j  *  1,10  -  angular  index 

N  >  1  -  number  of  onlecules  (spec  1) 

N  »  2  -  number  of  molecules  (spec  2) 

N  =  3  -  number  of  molecules  (spec  3) 

(starting  address  -  1)  of  molecules  as  ordered  in  (IP) 


N  =•  4 


Reg(N,lcR,kS)  Data  table  for  a  specific  region  (real) 


kR  -  1,10 
kS  -  1,20 
N  -  1 


N  -  2 
N  -  3 
N  -  4 
N  -  5 


index  of  a  region  in  a  sector 
index  of  Che  sector 

Dfl  »  differential  angle  (axisymmetric) 

DPI  is  a  weighting  factor 
DNl  «  number  density  (species  1) 

DN2  «  number  density  (species  1) 

VOL  1  “  actual  volume  of  smallest  cell  in  kR 
AREAl  IllPUT  area  of  smallest  cell  in  kR 


Region  Geometry  and  Input  Flux 


F1(10),F2(10) 


(equal  number  of  molecules  of  either  jet  species) 
input  flux  through  high  pressure  starting  line 
(spec  1),  (spec  2) 


Input  Flux  From  High  Pressure  Starting  Line 


FWPl(N,I,j) 
FWP2(N,I,j) 
I  I  I  species 


(positive  (input  molecules) 


I  “  1,10  -  number  of  the  cell  along  the  starting  line  in  a  region 


J  =  1,10  -  number  of  the  region  along  the  starting  line 


N  *  1  -  molecular  flow  for  a  given  DFl 

N  *  2  -  mean  molecular  velocity  Vx 

N  =  3  -  mean  molecular  velocity  Vy 

N  =»  4  -  mean  input  gas  temperature 


7 .  Output  Flux 

FOTn(4,NAD,kR) 

FNN2(4,NAD,kR) 

•1  + 

I  negative  (output) 

I 

north 


FSNl(4,NAD,kR) 

FSN2(4,NAD,kR) 

negative  (output) 
south 


kR  -  number  of  region  in  the  sector 
NAD  -  angular  location 

The  first  index  include  the  same  parameters  as  (FWPl) 


FEN1(4,NRD) 
FEN2(4,NRD) 
I  negative 


FWN1(4,NRD) 

FWN2(4,NRD) 

ft 

I  negative 
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NRD  radial  location 


FENl ,FEN2 ,FWN1 ,FWN2  are  necessary  for  Iterations  within  one  sector 
Sampling  of  Output  Flux  from  a  Sector 

After  averaging  they  are  transferred  to  FOEl ,  F0E2,  FOWl ,  F0W2. 
F0Wl(4,kC,kR,kS) 

output  flow  to  the  west 

F0W2(4,kC,kR,kS) 

F0El(4,kC,kR,kS) 

output  flow  to  the  east 

FOE2(4,kC,kR,kS) 

kC  -  radial  location  (cell)  of  flow  In  a  region 
kR  ~  Index  of  a  region  In  a  sector 
kS  -  Index  of  the  sector 

The  first  Index  (4  parameters)  Include  the  same  parameters  as  FEP1(  ) 


B.2  MOLECULAR  SIMULATION  FOR  A  GIVEN  REGION 

After  we  define  the  geometry  of  the  whole  sector  resulting  from  the 
region  geometry  and  cell  geometry,  we  may  start  with  the  molecular  simulation. 
This  Includes : 

-  Initial  setting  of  molecules  In  cells 

-  molecules  are  moved  according  to  the  time  Increment  DTl-1 

-  new  molecules  are  generated  according  to  Input  (or  output)  flows 

-  collisions  calculations 

-  Integration  of  flow  parameters  for  average  parameters  calculations 

-  repetition  of  the  whole  procedure  as  long  as  necessary  to  obtain 
reasonable  statistical  averaging 

-  calculate  averages  and  flow  weighting. 

These  routines  are  the  core  of  the  program  and  must  be  repeated  for  all 
regions  In  a  sector  and  (or  all  sectors  In  the  domain  where  the 
collisions  are  significant). 

In  the  following  sections  we  bring  a  detailed  description  of  this  part  of 
the  program. 

1 .  Initial  Setting  of  Molecules  In  Cells 

a.  The  Initial  number  of  available  molecules  In  a  simulation  is 
larger  than  the  number  of  active  molecules  (PI (data,  number  of 
molecules),  P2(data,  number  of  molecules)  are  the  vectors  used  for 
species  1  and  2) 

b.  an  Inactive  molecule  is  defined  as 


[P1(4,N)  or  P2(4,N)]  =  -99 


c.  calculation  of  number  of  molecules  to  be  set  in  each  cell 

d.  calculation  of  cell  coordinates 

e.  deactivation  of  all  available  molecules  In  simulation 

f.  definition  of  molecules  coordinates 
P1(4,N),  P2(4,M)  are  polar  radiuses 

P1(5,N),  P2(5,M)  are  angular  coordinate  in  radians 

All  molecules  in  a  cell  are  set  at  random  locations  within  the 
cell. 

g.  definition  of  molecular  velocities 
P1(1,N)»  P2(1,M)  velocity  in  X  direction 
P1(2,N),  P2(2,fl)  velocity  in  Y  direction 
P1(3,N),  P2(3,M)  velocity  in  Z  direction 

Thermal  velocities  are  random  function  of  temperatures  and  are 
added  to  the  mean  velocity  as  defined  at  initial  boundary  ALFAq  of 
the  region. 

As  the  thermal  velocity  has  a  Boltzmann  distribution  the  thermal  velocity 
setting  is  based  on  rejection-acceptance  methods  (for  more  details  see  Bird 
[4J  Appendix  D) . 


h.  reset  collision  timers  and  relative  velocity 


1.  reset  general  time  counter:  Time  =  0 
The  Simulation 

a.  Move  all  molecules  according  to  their  velocity  (V^,  Vy)  and  find 
their  new  coordinates 

Note  -  A  routine  designed  to  calculate  the  collisions  with  the 
wall  was  included  in  the  program;  if  the  region  (or  sector)  is 
bounded  by  the  solid  wall,  the  collision  is  calculated  -  resulting 
new  velocities  and  directions  and  counted  for  wall  flux 
calculations.  If  the  program  is  stopped  at  an  angle  where  the 
flow  becomes  collisionless,  this  routine  becomes  irrelevant  and 
other  type  of  calculations  should  be  designed. 

b.  Output  flow  counting:  all  molecules  that  leave  the  region  are 
counted  and  stored  in  specific  vectors  which  are  used  as 

inputs  to  other  regions.  The  output  vectors  are  (X  represent  1  or 
2  for  the  two  species  in  the  program  (FSNX)) 

F^NX(l,j,kR)  •>  "south"  boundary 
FN^NXC  1 , j  ,kR)  ->•  "north"  boundary 
F^NX(1,I)  "east"  boundary 

F^NX(1,I)  +  "west"  boundary 

I  represents  the  radial  location  index  of  the  cell 
j  represents  the  angular  location  index 
kR  is  the  index  of  the  region  within  a  sector 


Note  -  all  molecules  that  move  to 

(j=NAIH-l)(I-URD+l)  are  placed  in  FENX(1,NRD) 

(j=  -1)(X=NRIH-1)  are  placed  in  FWNX(1,NRD) 

(j=NAIH-l)(I-  -1)  are  placed  in  FSNX(l,l,kR) 

(j=  -1)(I=  -1)  are  placed  in  FWNX(l,l,kR) 

This  was  done  only  for  simplification  reasons. 

Generation  of  new  molecules  due  to  input  flows.  Through  the  four 
sides  E,N,W,S,  of  a  region,  molecules  are  allowed  to  enter  the 
region  according  with  the  flows  coming  from  the  neighbouring 
regions : 

♦ 

for  "W"  side  of  the  region  in  sector  1 

FWPX(P,I,kR) 

for  other  sectors 


for  "E"  side  of  any  region 
for  "N"  side  of  any  region 
for  "S"  side  of  the  region 


FOEX(P,I,kR,kS-l) 
FOWX(P,I,kR,kS+l) 
FSNX(P,j ,kR+l) 
FNNX(P,j ,kR-l) 


The  first  parameter  of  all  these  arrays  represent: 
P  =  1  -  number  flux  (real  number) 

P  =  2  -  velocity  component  -  (Vx) 

P  =  3  -  velocity  component  -  (Vy) 


P  =  4  -  gas  temperature 


Note  1  -  because  every  region  has  a  different  size  of  £ingle  DFI,  all 
fluxes  have  to  be  adjusted  accordingly. 

Note  2  -  Input  fluxes  are  calculated,  adjusted  and  stored  as  real 
numbers.  The  number  of  Input  molecules  are  by  definition  Integers. 

In  order  not  to  "loose"  molecules,  the  number  of  Input  molecules  Is 
Increased  by  1  on  a  random  basis.  (The  average  of  many  runs  will 
result  In  the  accurate  average  Input  flow.) 

New  molecules  are  set  at  random  locations  on  the  boundary  of  the  specific 
cells  and  at  random  time  within  DTli.  Then  each  molecule  is  allowed  to  enter 
the  region  according  to  its  Initial  coordinate  and  velocity.  At  the  end  of 
the  time  interval  the  new  location  and  velocity  is  stored  in  molecule  array 
PI  or  P2.  If  DTM  is  chosen  to  be  too  large  and  cell  size  is  small  (total 
region  size  too  small)  some  molecules  may  cross  the  region  and  will  not  be 
counted  in  the  simulation  of  the  specific  region.  In  order  not  to  "loose" 
molecules : 

(a)  DTM  should  be  decreased 

(b)  count  these  molecules  as  output  fluxes  from  the  specific 
region.  (This  is  recommended  only  if  there  is  no  other 
choice . ) 

Note  -  because  the  arrays  of  input  flows  store  only  averagea  data  for 
the  molecules ,  the  thermal  velocity  of  each  new  molecule  is  calculated 
according  to  the  Boltzmann  distribution  as  a  function  of  the  averaged 
temperature . 

d.  Rearrangement  of  molecules  in  cells.  Before  collisions  are 

calculated  all  simulated  molecules  which  have  been  let  to  raove  and 
generated  have  to  be  rearranged  and  recounted  for  each  cell. 


The  array  IC(Ic,I,J)  contains  Integer  data  for  each  cell 

(I.j) 

k  1  -  is  the  number  of  molecules  of  species  1 

k  ■  2  -  is  the  number  of  molecules  of  species  2 

k  ■  3  -  is  the  number  of  molecules  of  species  3  (not  used) 

k  ■  4  -  is  the  (address-1)  of  the  first  molecule  in  the  cell 

related  to  the  vector  IP(M) 

Vector  IP(M)  contains  the  list  of  the  simulated  molecules  arranged 
in  the  order  of  species  in  cells  and  cells  respectively.  The 
following  is  a  graphic  description  of  IP(M) 

rS'A/  rs 


Figure  23.  Vector  IP(Il) 


B.3  SUlUL  Program  Flowchart 

A  simplified  flowchart  for  the  Monte  Carlo  simulation  of  the  molecular 
flow  is  shown  in  Figure  24.  The  program  is  designed  to  solve  the  ring 
axlsymmetrlc  jet  flow,  however,  minor  changes  may  be  done  to  enable  a 
different  geometry. 


B.4  SIMUL  Program  Listing 

C  PROGRAM  SIMUL  SIMOOOlO 

C  THIS  PROGRAM  IS  DESIGNED  TO  CALCULATE  THE  MOLECULAR  FLOW  OUTSIDE  THE  SIM00020 

C  CONTINUUM  REGION  FLOW  OF  A  HIGHLY  UNDEREXPANDED  AXISYMMETRIC  RING  JET.SIM00030 
C  RESULTS  FOR  THE  CONTINUUM  FLOW  MAY  BE  OBTAINED  FROM  'AXSYM'  PROGRAM  SIM00040 

C  WHICH  GIVES  THE  ‘METHOD  OF  CHARACTERISTICS'  ISENTROPIC  SOLUTION.  SIM00050 

C  THE  BOUNDARY  BETWEEN  THE  CONTINUUM  AND  MOLECULAR  FLOW  IS  DEFINED  BY  SIM00060 

C  THE  BREAKDOWN  PARAMETER  'P'  AS  PROPOSED  BY  G.  A.  BIRD.  SIM00070 

C  THE  'MOLECULAR  DOMAIN'  IS  DIVIDED  INTO  POLAR  DIVISIONS  MAKING  A  SET  SIMQOOSO 

C  OF  SECTORS.  EACH  SECTOR  IS  SUDIVIDED  INTO  10  REGIONS  SIM00090 

C  NO  APRIORY  INFORMATION  ABOUT  THE  GEOMETRY  OF  THE  DIFFERENT  SIMOOlOO 

C  REGIONS  IS  AVAILABLE,  THEREFORE,  MANUAL  INTERVENTION  MAY  BE  SIMOOllO 

C  REQUIRED  WHEN  MOVING  FROM  ONE  SECTOR  TO  THE  OTHER.  AN  ACCEPTABLE  SIM00I20 

C  GEOMETRY  WILL  RESULT  A  REASONABLE  NUMBER  OF  SIMULATED  MOLECULES.  SIM00130 

C  EACH  REGION  IS  SUBDIVIDED  INTO  A  NUMBER  OF  CELLS  WITH  A  GEOMETRY  SIM0D140 

C  DEFINED  BY  A  POLAR  MESH.  A  NUMBER  OF  MOLECULES  IS  SET  IN  EACH  CELL  SIM00150 

C  PROPORTIONAL  TO  THE  CELL  VOLUME.  THE  BOUNDARY  CONDITIONS  FOR  EACH  SIM00160 

C  REGION  REQUIRE  INPUT  AND  OUTPUT  FLUX  OF  MOLECULES.  NO  APRIORY  SIM00170 

C  INFORMATION  ON  THE  FLUX  IS  AVAILABLE.  IT  WILL  BE  CALCULATED  IN  SIM00180 

C  AN  ITERATIVE  MODE.  SIM00190 

C  IF  THE  BOUNDARY  CONDITIONS  FOR  ALL  CELLS  ARE  CONSTANT  THE  NUMBER  SIM00200 

C  FLUX  INTO  EACH  CELL  IS  PROPORTIONAL  TO  CELL  WALL  AREA.  SIM00210 

C  TO  DECREASE  THE  ERROR  WHEN  INTRODUCING  MOLECULES  -(INTEGER  NUMBER)-  SIM00220 

C  DUE  TO  THE  INPUT  FLUX  -(REAL  NUMBER),  AN  ADDITIONAL  MOLECULE  IS  SIM00230 

C  GENERATED  ON  THE  BASIS  OF  RANDOM  NUMBERS  SUCH  THAT  THE  AVERAGE  OF  A  SIM00240 

C  LARGE  NUMBER  OF  SAMPLINGS  WILL  EQUAL  THE  REAL  INPUT  FLUX.  SIM00250 

CXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  SIM00260 


THE  MONTE  CARLO  SIMULATION 
THE  JET  GAS  IS  COMPOSED  OF  TWO  SPECIES  OF  MOLECULES 
AMBIENT  GAS  IS  REGARDED  AS  ONE  SPECIES 
THE  MOLECULAR  MODEL  IS  -  'HARD  SPHERE  MOLECULE' 
NETWORK  DEFINITION 


SIM00270 

SIM00280 

SIM00290 

SIM00300 

SIM003I0 


THE  MAXIMUM  RADIUS  (POLAR)  OF  THE  DOMAIN  IS  ASSUMED  TO  BE  RP=15  M.  SIM00320 
THE  ANGLE  OF  THE  BOUNDARY  BETWEEN  CONTINUUM  AND  MOLECULAR  FLOW  SIM00330 
AND  THE  SOLID  WALL  IS  SIM00340 

ALFA  (CALCULATED  IN  PROGRAM  'AXSYM')  SIM00350 

DEFINE  A  POLAR  SECTOR  WITH  A  RADIUS  'RP'  AND  AN  ANGLE  OF  SIM00360 

DALFA  =5XMFPXNAD/RP  SIM00370 

THIS  SECTOR  IS  SUBDIVIDED  INTO  A  NUMBER  OF  RADIAL  DIVISIONS  SIM00380 

MAKING  'N'  REGIONS  FOR  SIMULATION  CALCULATIONS  FOR  EACH  'DALFA'.  SIM00390 
EACH  REGION  IS  DEVIDED  INTO  SIM00400 

NAD  -ANGULAR  DIVISIONS  (15)  WITH  AN  ANGLE  OF  DDALFA=5XMFP/RP  SIM00410 

NRD  -RADIAL  DIVISIONS  (10)  SIM00420 

MAKING  NADXNRD  CELLS  SIM00430 

THE  SMALLEST  CELL  CONTAINS  'MIN'  MOLECULES  SIM00440 

MIN  =  15  SIM00450 

TOTAL  NUMBER  OF  ACTIVE  MOLECULES  IN  A  REGION  IS  LIMITED  TO  6000  SIM00460 

(3000  MOLECULES  OF  EACH  SPECIES.)  SIM00470 

THE  WIDTH  OF  A  CELL  DEFINED  BY  ANGLE  'DFI'  MAY  NOW  BE  EVALUATED.  SIM00480 

(MIN/NUMBER  DENSITY  =  ACTUAL  CELL  VOLUME)  SIM00490 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00500 
DEFINE  PARAMETERS . BLOCK  1SIM00510 


COMMON  IX 

DIMENSION  P1(5,3000),P2(5,3000),IP(6000),C(20,10,15), 
XIC(4.10,15),REG(5, 10,20) 

DIMENSION  DFI(10),DN(10) 

DIMENSION  R(10),A(10),VOL(10),M1(10),F1(10),F2(10) 

DIMENSION  SS(2,10,10) 

DIMENSION  FEN1(4,10),FWN1(4,10),FNN1(4,15,10),FSN1(4,15,10) 
DIMENSION  FEN2(4,10),FWN2(4,10),FNN2(4,15,10),FSN2(4,15,10) 
DIMENSION  FWP1(4,10,10),FWP2(4,10,10) 

DIMENSION  FOE1(4,10,10,20),FOW1(4,10,10,20) 

DIMENSION  FOE2(4,10,10,20),FOW2(4,10,10,20) 

DIMENSION  NM(3),VRC(3),TOC(3,3),SPEC(3,5).NCOL(10,15) 

P1,P2,P3  CONTAIN  INFORMATION  ON  UP  TO  M  SIMULATED  MOLECULES 

P1(1,N),P1(2,N),P1(3.N)  ARE  U,V,W  VELOCITY  COMPONENTS  (CARTESIAN) 
P1(4,N),P1(5,N)  ARE  R  AND  TETA  COORDINATES  (POLAR) 


SIM00520 

SIM00530 

SIM00540 

SIM00550 

SIM00560 

SIM00570 

SIM00580 

SIM00590 

SIM00600 

SIM00610 

SIM00620 

SIM00630 

SIM00640 

SIM00650 

SIM00660 

3IM00670 

SIM00680 

SIM00690 


IP(M)  ARE  THE  M  MOLECULES  ARRANGED  IN  THE  ORDER  OF  THEIR  CELLS 
C  C(20,I.J)  CONTAINS  INFORMATION  ON  UP  TO  Ixj  CELLS 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXKXKXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXSIM00700 
C  SIM00710 

C  ALFA  IS  THE  ANGLE  WHERE  THE  BREAKDOWN  PARAMETER  EQUALS  .05  SIM00720 


78 


oooo 


C  FPM  IS  THE  MEAN  FREE  PATH  AT  ALFA.  SIM00730 
CXXXXXXXKXXXXXXXXXXXXXXXXifXKXXXXXJfIfXXKXXXXXXXXXXXXXXXliXXXXXXXXXXXXXXXXXXSIMOOTAO 
C  SET  GENERAL  CONSTANTS . BLOCK  2SIM00750 


IX 
PI 

BOLTZ 

PO 

TO 

AV06 
R6 


529814367 
3.141593 
1.38044E-23 
101325. 

273. 

2.68699E^25 
8314. 


SIM00760 

SIM00770 

SIM00780 

SIM00790 

SIM00800 

SIM00810 

SIM00820 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM00830 

C  SET  PROGRAM  CONSTANTS . BLOCK  3SIM00840 

ISPEC=2  SIM00850 

R1  -  2.5  SIM00860 

DR3.15  SIM00870 

RP-15.  SIM00880 

C  NUMBER  OF  SIMULATED  MOLECULES  IN  SMALLEST  CELL  (DIFFERENT  SPECIES)  S1M00890 

NMOLl-3000  SIM00900 

NMOL2=3000  SIM00910 

NMOL3=0  SIM00920 

MIN  =  15  SIM00930 

NIM  -  15  SIM00940 

N2M  =  15  SIM00950 

C  NUMBER  OF  DIVISIONS  IN  A  SIMULATED  REGION  SIM00960 

NAD  »  15  SIM00970 

NRDS^IO  SIM00980 

NRD=10  SIM00990 

NRDl-  3  .  SIMOIOOO 

NIS  =2  SIMOlOlO 

C  R1  IS  THE  RADIUS  OF  THE  CYLINDER  (MALL)  SIM01020 

C  R(I)  IS  THE  RADIAL  COORDINATE  OF  CELL(I)  SIM01030 

C  VOL(I)  IS  THE  RATIO  BETWEEN  VOLUMES  OF  CELL(I)  AND  CELL(l)  SIM01040 

C  A(I)  IS  THE  RATIO  BETWEEN  INPUT  FLOW  AREAS  OF  CELL(I)  AND  CELL(l)  SIM01050 

C  DR  IS  THE  RADIAL  SIZE  OF  A  CELL  (RADIAL  INCREMENT)  SIM01060 

C  INPUT  THE  FOLLOWING  AS  ‘DATA*  OR  'READ*  STATEMENTS  xxxxxxxxxxxxxxxx  SIM01070 

TETA  =1.2  SIM01080 

C  TETA  IS  THE  FLOW  DIRECTION  (RADIANS)  ON  THE  BREAKDOWN  LINE  AS  FOUND  SIM01090 

C  FROM  THE  AXSYM  PROGRAM.  VO  IS  THE  FLOW  VELOCITY  (M/SEC)  SIMOllOO 

VO  =  2100  SIMOlllO 

VOX  =  VO»COS(TETA)  SIM01120 

VOY  =  VO*SIN<TETA)  SIM01130 

TWALL=SOO.  SIM0I14D 

DN01=1.1E21  SIM01150 

DN02=1.1E21  SIM01160 

DTM=1.5E-6  SIM01170 

C  SET  DATA  FOR  THE  DIFFERENT  SPECIES  SIM01180 

C  MOLECULAR  MASS  SIM01190 

SPEC(1,1)=4./AV0G  SIM01200 

SPEC(2,1)=40./AV0G  SIM01210 

SPEC(3,1)=29./AV0G  SIM01220 

C  MOLECULAR  DIAMETER  SIM01230 

SPEC(1,2)=2.19E-10  SIM01240 

SPEC(2,2)=4.00E-10  SIM01250 

SPEC(3,2)=  SIM01260 

SPEC(1,3)=  SIM01270 

SPEC(2,3)=  ADDITIONAL  DATA  IF  REQUIRED  SIM01280 

SPEC(3,3)=  SIM01290 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01300 


C  THERMAL  VELOCITIES  AT  WALL  TEMPERATURE 
VWM1=SORT(2.XBOLTZXTWALL/SPEC(1,1)) 
VWM2=S0RT(2.»B0LTZ»TWALLXSPEC(2,1)) 
C 


SIM01310 

SIM01320 

SIM01330 

SIM01340 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXKKXXXKXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01350 


C  INITIALISATION. .RESET  ALL  SAMPLING  VARIABLES 
DO  80  1=1, NRD 
DO  80  JR=1,10 
DO  80  JS=1,20 
DO  80  KPAR=1,4 
F0E1(KPAR,I,JR,JS)=0. 

F0E2(KPAR,I,JR, JS)=0. 

F0W1(KPAR,I,JR, JS)=0. 

80  F0W2(KPAR,1, JR, JS)=0. 


SIM01360 

SIM01370 

SIM01S80 

SIM01390 

SIM01400 

SIM01410 

SIM01420 

SIM01430 

SIM01440 


CXXXXXXXlIXXXXXXXXafXXXXXXXXXXXXXXKXlIKlIXlIXKXXKXXXXXXXXXXXXXXXXKXlfXXifXXXBfXXSIMOl^SO 

C  SET  INPUT  PARAMETERS . BLOCK  S  SIM01460 

ITER=1  SIM01470 

C  RETURN  TO  3000  FOR  ADDITIONAL  ITERATIONS  SIM01480 

3000  KS-1  SIM01490 

ALFA*1.3  SIM01500 

FPM=.001  SIM01510 

TEMPS40.  SIM01520 

VTER1=SQRT(2. XBOLTZXTEMP/SPECCl . 1 ) )  SIM01530 

VTER2-SQRT( 2 . XB0LTZXTEMP/SPEC(2, 1 ) )  SIM01540 

IF  (ITER.GT.DGO  TO  2005  SIM01550 

DO  2001  1-1,10  SIM01560 

REG(2,I,1)=DN01  SIM01570 

REG(3,I,1)=DN02  SIM01580 

DO  2001  J-1,10  SIM01590 

FNP1(2,I,J)-V0X  SIM01600 

FWP2(2,I,J)=V0X  SIM01610 

FWP1(3,I,J)=V0Y  SIM01620 

FWP2(3,I,J)=V0Y  SIM01630 

FWP1(4.1,J)=TEMP  SIM01640 

2001  FMP2(4,I,J)=TEMP  S1M01650 

2005  CONTINUE  SIM01660 

C  SIM01670 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX9fXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01680 
C  SIM01690 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01700 
C  DEFINE  SECTOR  GEOMETRY . BLOCK  4  SIM01710 


C  KR  IS  THE  INDEX  FOR  A  REGION  IN  ONE  SECTOR  SIM01720 

C  RETURN  TO  1000  FOR  NEXT  SECTOR  SIM01730 

DALFA=0.  SIM01740 

1000  KR=1  SIM01750 

C  RESET  SAMPLING  OF  FLOM  VARIABLES  (N  AND  S)  SIM01760 

DO  81  1=1, NAD  SIM01770 

00  81  JR=1,10  SIM01780 

DO  81  KPAR=1,4  SIM01790 

FNNKKPAR,!,  JR)  =  0.  SIM01800 

FNN2(KPAR,I,JR)=0.  SIM01810 

FSNl(KPAR,I,JR)sO.  SIM01820 

81  FSN2(KPAR,I,JR)=0.  S1M01830 

C  ALFA  IS  THE  STARTING  ANGLE  OF  THE  PRESENT  SECTOR.  PREVIOUS  DALFA  SIM01840 

C  SUBTRACTED  FROM  PREVIOUS  ALFA  (PREVIOUS  SECTOR)  GIVES  THE  NEW  ALFA..  SIM01850 

ALFA=ALFA-DALFA  SIM01860 

DALFA  =  5.XFPMXFL0AT(NAD)/RP  SIM01870 

C  DALFA  IS  THE  ANGLE  OF  THE  SECTOR  SIM01880 

IF(DALFA. GT. ALFA/2. )DALFA=ALFA/2.  SIM01890 

IF( DALFA. GT.(FPMX75/RP))DALFA=ALFA  SIM01900 

IFCALFA.LE.O. )G0  TO  2000  SIM01910 

DDALFA  =  DALFA/FLOAKNAD)  SIM01920 

C  DDALFA  IS  THE  ANGLE  OF  A  SINGLE  CELL  IN  THE  SECTOR  SIM01930 

ALFAJ=ALFA-DALFA/2.  SIM01940 

C  SIM01950 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM01960 

C  DEFINE  REGION  KR  IN  SECTOR . BLOCK  5  SIM01970 

C  RETURN  TO  2000  FOR  THE  NEXT  REGION  KR  SIM01980 

2000  CONTINUE  SIM01990 

C  RESET  SAMPLING  OF  FLOW  VARIABLES  PER  REGION  SIM02000 

DO  82  1=1, NRO  SIM02010 

DO  82  KPAR=1,4  SIM02020 

FEN1(KPAR,I)=0.  SIM02030 

FEN2(KPAR,I)=0.  SIM02040 

FWN1CKPAR,I)=0.  SIM02050 

82  FWN2(KPAR,I)=0.  SIM02060 

MT=0  SIM02070 

NRD=NRDS  SIM02080 

IF(KR.LT.2)NR0=NRD1  SIM02090 

DO  100  I  =  1,NRD  SIM02100 

C  POLAR  RADIUS  MEASURED  FROM  THE  NOZZLE  LIP  SIM02110 

R(I)  =(FL0AT(I)-.5)XDR  SIM02120 

IF(KR.EQ.2)R(I)=R(I)+FL0AT(NRD1)XDR  SIM02130 

IF(KR.GE.3)R(I)=R(I)+(FL0ATCNRD1+(KR-2)XNRD))*DR  SIM02140 

C  RSI  IS  THE  RADIUS  MEASURED  FROM  THE  AXIS  OF  SYMMETRY  SIM02150 

RSI=R(I)+R1/SINCALFAJ)  SIM02160 


c. 

c 


IF(I.EQ.1)RS1=RSI 

A(I)  -  RSI/RSl 

VOL(I)  »  R(I)XRSI/(R(1)XRS1) 

Ml(I)  s  MINXVOL(I) 

MT  =  MT+M1(I) 

NIO  IS  THE  INITIAL  NUMBER  OF  MOLECULES  IN  EACH  CELL 
DN1-REG(2,KR,KS) 

DN2=REG(3,KR,KS) 

REG(1,KR,KS)SFL0AT(MIN)/(DN1XR(1)XRS1XDRXDDALFA) 

DFI  (REG(1,KR,KS))  IS  THE  SPHERICAL  (AXISYMMETRIC)  ANGLE 
THIS  ANGLE  HAS  DIFFERENT  VALUES  FOR  EACH  KR,  THEREFORE  IT  IS  A 
WEIGHTING  FACTOR . +=+=+=+=+*+*+8+=+=+=+=+*+=+=+=+=+=+=+*+++++++++++++ 
(FNNI,FNN2,FSN1,FSN2.FEN1,FEN2,FWN1,FWN2)  X  DFI 

(FOE1,FOE2,FOW1,FOW2)  X  DFI 

DAI  =  RSlXDRXREGd.KR.KS) 

DAI  IS  THE  INPUT  AREA  OF  CELL(l) 

FNl  =  V0XDA1XREG(2,KR,KS)XSIN<ALFA-TETA) 

FN2  s  VOXOA1XREG(3.KR,KS)XSIN(ALFA-TETA) 

TETA  IS  THE  ANGLE  BETWEEN  FLOW  DIRECTION  AND  THE  WALL 
FKI)  =  FNIXA(I) 

F2(I)  =  FN2XA(I) 

99  FWP1(1,I,KR)=F1(I) 

100  FWP2(1,I,KR)=F2(I) 

DEHNE  ‘  ACTUAL  ‘  VOLUME  ’  AND 'input ‘area ‘6f‘ SMALLEST ‘cell 'in  REGION 
REG(4.KR,KS)-R(1)XDDALFAXDRXRS1XREG(1,KR,KS) 

REG(5,KR,KS)-DA1 

FNl  IS  THE  NUMBER  FLUX  TO  THE  SMALLEST  CELL  (REAL  NUMBER)  PER  SECOND. 
INPUT  NUMBER  OF  MOLECULES  (INTEGER)  WILL  BE  INTEGRATED  TO  MAKE  AN 
AVERAGE  OF  FNl. 

FKI)  IS  THE  INPUT  FLUX  TO  CELL(I) 

FOLLOWING  ARE  THE  REGION  BOUNDARIES 

102  RMIN-R(l)-.SxDR 
RMAX=RMIN-t-DRXFLOAT(NRD) 

TMAX=ALFA 

TMIN=ALFA-DALFA 

DO  9102  I  =  1,NRD 

WRITE  (6,9101)R(I),A(I),VOL(I),M1(1),F1(I) 

FORMATC  SSFlS.SdlO^ElS.S) 

CONTINUE 
MTT=  MTXNAD 
WRITE(6,103)MT,MTT 

103  FORMAT  (•  INITIAL  NUMBER  OF  SIMULATED  MOLECULES  IS=  M5, '  PER  SP 
XCIES  PER  DDALFA,  TOTAL  NUMBER  IS',  15) 


9101 

9102 


C 

C 

c 


SIM02680 
SIM02690 
SIM02700 
SIM02710 
SIM02720 
SIM02730 
SIM02740 
SIM02750 
SIM02760 
SIM02770 
SIM02780 
SIM02790 
SIM02800 
SIM02810 
SIM02820 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXJOfXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXX  SIM02830 

SIM02840 
SIM02850 
SIM02860 
SIM02870 
SIM02880 


C 

C 


CALCULATE  INITIAL  NUMBER  OF  MOLECULES  IN  CELLS . BLOCK  7 

SET  INITIAL  STATE  OF  GAS .  "  8 

DO  150  1  =l.NRO 
DO  150  J=  l,NAO 

SET  SIMULATED  MOLECULES  IN  THEIR  CELLS 
IC(1,I,J)  =  MKI) 

IC(2,I,J)  =  MKI) 

IC(3,I,J)  =  0 

SET  CELL  COORDINATES 
C(19,I,J)  =  R(I) 

C(20,I,J)  =  ALFA-(FL0AT(J)-.5)XDDALFA 
150  CONTINUE 


C  DEACTIVATE  ALL  MOLECULES 
DO  170  N  -  1,NM0L1 
PK4,N)  s  -99. 

170  P2(4,N)  »  -99. 


81 


or>  o  oo  o  o  DO 


C  SIN0Z890 

CXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXKKXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXafXX  SIM02900 
SET  INITIAL  STATE  OF  THE  GAS  (LOCATION  AND  VELOCITY  OF  MOLECULES)  SIM02910 

SIM02920 


NADRl  »  0  SIM02930 

NADR2  -  0  SIM02940 

DO  200  I  =  1,NRD  SIM02950 

DO  200  J  »  I, NAD  SIM02960 

NMl  =  IC(1,I,J)  SIM02970 

SIM02980 

DO  205  N  =  1,NM1  SIM02990 

NADRl  -  NADRl  +  1  SIM03000 

CALL  RANDU(PP)  SIM03010 

PKA, NADRl)  =  C(19,I,J)+DRX(PP-.5)  SIM03020 

CALL  RANDU(P)  SIM03030 

Pl(5, NADRl)  =  C(20.1,J)-t-DDALFAX(p-.5)  SIM03040 

SIM03050 

DO  205  NV  =  1,3  SIM03060 

203  CALL  RANDU(P)  SIM03070 

V  =  -3.+6.XP  SIM03080 

B  =  EXP(-VXV)  SIM03090 

CALL  RANDU(P)  SIM03100 

IF(B.LT.P)  GO  TO  203  SIM03110 

PKNV. NADRl)  =  PXSIN(B)XVTER1  SIM03120 

IF(NV.E0.1)P1(NV,NADR1)=P1(NV,NADR1)+V0X  SIM03130 

IF(NV.EQ.2)P1(NV,NADR1)=P1(NV,NADR1)+V0Y  SIM03140 

205  CONTINUE  SIM03150 

SIM03160 

REPEAT  PROCEDURE  FOR  SPECIES  2  SIM03170 

NM2  =IC(2,I,J)  SIM03180 

DO  210  N  =  1,NM2  SIM03190 

NADR2  =  NADR2-t-l  SIM03200 

CALL  RANDU(P)  SIM03210 

P2(4,NADR2)  =  C(19,I, J)+DRx(P-.5)  SIM03220 

CALL  RANDU(P)  S1M03230 

P2(5,NADR2)  =  C(20,I, J)+DDALFAX(P-.5)  SIM03240 

SIM03250 

DO  210  NV  *  1,3  SIM03260 

207  CALL  RANDU(P)  SIM03270 

V  =  -3.+6.XP  SIM03280 

B  =  EXP(-VXV)  SIM03290 

CALL  RANDUCP)  SIM05300 

IF(B.LT.P)  GO  TO  207  SIM03310 

P2(NV,NADR2)  =  PXSIN(B)XVTER2  SIM03520 

IF(NV.EQ.1)P2(NV,NADR2)=P2<NV,NADR2)+V0X  SIM03330 

IF<NV,EQ.2)P2(NV,NADR2)=P2(NV,NADR2)+V0Y  SIM03340 

210  CONTINUE  SIM03350 

SIM03360 

WHEN  NECESSARY,  REPEAT  PROCEDURE  FOR  SPECIES  3.  SIM03370 

200  CONTINUE  SIM03380 

C  SIM03390 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03400 


CALL  RANDU(P) 
P2(5,NADR2)  = 


C(20,I,J)+DDALFAX(p-.5) 


DO  210  NV  *  1,3 
207  CALL  RANDU(P) 

V  =  -3.+6.XP 
B  =  EXP(-VXV) 

CALL  RANDU(P) 

IF(B.LT.P)  GO  TO  207 
P2(NV,NADR2)  =  PXSIN(B)XVTER2 
IF(NV.EQ.1)P2(NV,NADR2)=P2<NV,NADR2)+V0X 
IF<NV,EQ.2)P2(NV,NADR2)=P2(NV,NADR2)+V0Y 
210  CONTINUE 

WHEN  NECESSARY,  REPEAT  PROCEDURE  FOR  SPECIES  3. 
200  CONTINUE 


C  DEFINE  HERE  ALL  COLLISION  PARAMETERS  TO  BE  INCLUDED  IN  SIMULATION 


SIM03410 


C  RESET  COLLISION  TIMERS 
DO  19  1=1, NRD 
DO  19  J=1,NAD 
DO  19  L=l,3 
DO  19M=1,S 
KT=3X(L-l)+M+9 
C(KT,I,J)=0. 

C  SET  EXPECTED  MAXIMUM  RELATIVE  VELOCITY  IN  COLLISIONS 
KV=3X(L-1)+M 
EM1=SPEC(L,1) 

EM2=SPEC(M.l) 

EMR = EMI X  EM2/ ( EMI + EM2 ) 

C 

IF(KS.E0.1)G0  TO  17 
TEM=F0E1(4,I,KR,ICS) 

GO  TO  18 

17  TEM=FWP1(4,I,KR) 


BLOCK  9  SIM03420 
SIM03430 
SIM03440 
SIM03450 
SIM03460 
SIM03470 
SIM03480 
SIM03490 
SIM03500 
SIM03510 
SIM03520 
SIM03530 
SIM03540 
SIM03550 
SIM03560 
SIM03570 
SIM03580 
SIM03590 
S1M03600 


18  RV*2./SQRT(PIX2.XB0LTZ*TEM/EMR) 

19  C(KV,I,J)=2.KRV 


Siri03610 
SIM03620 
SIM03630 

C  MAXIMUM  RELATIVE  VELOCITY  HILL  BE  RESET  IF  FASTER  ENCOUNTERS  OCCUR  SIM03640 
C  SIM03650 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM03660 


TIME=0. 


SIM03670 


C  LOOP  OVER  TIME  INTERVALS . BLOCK  10  SIM03680 


C 

C 


DO  6000  JOTM  -  1,NIS 


TIME*TIME+DTM 


SIM03690 

SIM03700 

SIM03710 

SIM03720 


CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIN03730 
C  SIM03740 

C  MOVE  ALL  MOLECULES  . BLOCK  12  SIM03750 

C  MOVE  MOLECULES  OF  SPECIES  1  SIM03760 

C  SIM03770 

DO  310  II  -  l.NMOLl  SIMQ3780 

C  SKIP  INACTIVE  MOLECULES . 13  SIM03790 

IF  (P1(4,I1) .EQ.-99.)  GO  TO  310  SIM03800 

VX  =  PKl.Il)  SIM03810 

VY  =  P1(2,I1)  SIM03820 

RX  =  P1(4,I1)  SIM03830 

T  =  PKS.Il)  SIM03840 

TOLD=T  SIM03850 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIMQ3860 


C  FIND  NEW  COORDINATES . BLOCK  14 

XX=RXXCOS(T)+VXXDTM 
YY=RXXSIN(T)+VYXDTM 
RNEW=SQRT ( XXXX2+YYXX2 ) 

T=ATAN{YY/XX) 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  FOR  LAST  SECTOR  FIND  COLLISIONS  WITH  THE  WALL  AND  SAMPLE  THEM 

C  . BLOCK  15,16 

IF(ALFA.GT.DALFA)GO  TO  301 
IFCT.GT.O. )GO  TO  301 
DTR=DTMXT/<T-TOLD) 

C  DTR  IS  THE  TIME  REMAINING  AFTER  A  MOLECULE  STRIKES  THE  WALL 
IF(DTR.LT.1E-10)DTR=1E-10 
RW»RX+<RNEW-RX)XDTR/DTM 
IF(RW.LT.RMIN)RW=RMIN+DTRX.001 
I F  <  RW . GT . RMAX ) RW=RMAX- DTRX . 0  0 1 
L0C=RW/DR+1. 

C 

C  COUNT  COLLISIONS  WITH  THE  WALLCMUST  BE  WAIGHTED  =XDFI) 
SS(1,L0C,JSAMP)=SS(1,L0C,JSAMP)+1 

C 

C  SET  VELOCITY  AND  LOCATION  AFTER  A  MOLECULE  STRIKES  THE  WALL 
302  CALL  RANDU(P) 

B=VWM1XSQRT(-AL0G(P)) 

CALL  RANDU(P) 

BB=2.XPIXP 

P1(1,NADR1)=DXC0S(BB) 

VX=BXCOS(BB) 

P1(2,NADR1)=BXSIN(BB) 

VY=BXSIN(BB) 

CALL  RANDU(P) 

P1(3,NADR1)=VWM1XSQRT(-AL0G(P)) 

XX=RXXCOSCT)+VXXDTR 

YY=RXXSIN(T)+VYXDTR 

RNEW=S0RTCXXXX2+YYXX2) 

T=ATAN(YY/XX) 


SIM03870 
SIM03880 
SIN03890 
SIM03900 
SIM03910 
SIM03920 
SIM03930 
SIM03940 
SIM03950 
SIM03960 
SIM03970 
SIM03980 
SIM03990 
SIM04000 
SIM04010 
SIM04020 
SIM04030 
SIM04040 
SIM04050 
SIM04060 
SIM04070 
SIM04080 
SIM04090 
SIM04100 
SIM04110 
SIM04120 
SIM04130 
SIM04140 
SIM04150 
SIM04160 
SIM04170 
SIM04180 
SIM04190 
SIM04200 
SIM04210 
SIM04220 
SIM04230 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM04240 

C  DEACTIVATE  ALL  MOLECULES  THAT  MOVED  OUT . BLOCK  17  SIM04250 

.  301  IF(RNEW.LT. RMAX. AND. RNEW.GT.RMIN. AND. T.LT.TMAX. AND. T.GT.TMINIGO  T0SIM04260 

X303  SIM04270 

RI=(RNEW-RMIN)/DR+. 999999  SIM04280 

IR=RI  3IM04290 

TI=(TMAX-T)/DDALFA+. 999999  SIM04300 

IT=TI  SIM04310 

C  SIM04320 
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IF(IT.LE.O)GO  TO  304 
IF(IT.6T.NAD)60  TO  305 
IF(IR.6T.NR0)G0  TO  306 

C  COUNT  S  DIRECTION,  SAMPLE  . 

FSN1(1,IT,KR)  =  FSN1(1,IT,ICR)+1. 

GO  TO  309 

C  COUNT  M  DIRECTION,  SAMPLE . 

304  IF(IR.LE.0)IR=1 
IF(IR.GT.NRO)IR=NRD 
FNNl(l,IR)  =  FMNl(l,IR)-fI. 

GO  TO  309 

C  COUNT  E  DIRECTION,  SAMPLE . 

305  IF(IR.LE.0)IRsl 
IF(IR,GT.NRD)IR=NRD 
FEN1(1,IR)=FEN1C1,IR)+1. 

GO  TO  309 

C  COUNT  N  DIRECTION,  SAMPLE . 

306  FNNl(l,IT,KR)  =  FNNl(l,IT,KR)-l-l. 

SET  NEW  VALUES  IN  THE  MOLECULE  TABLE 


IT*  IT*  If"  'r*;' 


.18 

.18 

.18 

.18 


309  RNEW=-99. 

303  P1(4,I1)=RNEW 
P1(5,I1)=T 

310  CONTINUE 


REPEAT  PROCEDURE  FOR  SPECIES  2  MOLECULES. 

DO  320  I2=1,NM0L2 
SKIP  INACTIVE  MOLECULES 

IFCP2(4,I2) .EQ.-99. )GO  TO  320 
VX  =  P2(1,I2) 

VY  =  P2(2,I2) 

RX  =  P2(4.I2) 

T  =  P2(5,I2) 

T0LD=T 

XX=RXXC0S(T)+VX*DTM 

YY=RXXSIN(T)+VYKDTM 

RNEW=S0RT(XXX*2+YY»»2) 

T=ATAN(YY/XX) 

COLLISIONS  WITH  THE  WALL 


BLOCKS  12  T018 


IF(ALFA.GT.DALFA)GO  TO  311 
IFCT.GT.O. )G0  TO  311 
DTR=DTMKT/(T-TOLD) 

DTR  IS  THE  TIME  REMAINING  AFTER  A  MOLECULE  STRIKES  THE  WALL 
IF(DTR.LT.1E-10)DTR=IE-10 
RW=RX+ ( RNEW-RX ) XDTR/DTM 
IF(RW.LT.RMIN)RW=RW+DTRX.001 
I F ( RW . GT . RMAX ) RW= RW-DTRX .001 
LOC=RW/DR+l 

COUNT  COLLISIONS  WITH  THE  WALL  (MUST  BE  WEIGHTED  =XDFI) 

SS(2,L0C, JSAMP)=SS<2,L0C, JSAMP)+1 

FIND  THE  NEW  COORDINATES  OF  THE  MOLECULE 

312  CALL  RANDU(P) 

B=VWM2XS0RT(-AL0G(P)) 

CALL  RANDU(P) 

BB=2.xpixP 

P2(1,NADR2)=BXC0S(BB) 

VX=BXC0S(BB) 

P2(2,NADR2)=BXSIN(BB) 

VY=BXSIN(BB) 

CALL  RANDU(P) 

XX=RXXCOS(T)+VXXDTR 

YY=RXXSIN(T)+VYXDTR 

RNEW=S0RT(XXXX2+YYXX2) 

T=ATAN(YY/XX) 

DEACTIVATE  MOLECULES  THAT  MOVED  OUT. COUNT  FOR  OUTPUT  FLUX  EVALUATION 

311  IF(RNEW.LT. RMAX. AND. RNEW.OT.RMIN. AND. T.LT.TMAX. AND. T.GT.TMIN)GO 
X313 


SIM04330 
SIM04340 
SIM04350 
SIM04360 
SIM04370 
SIM04380 
SIM04390 
SIM04400 
SIM04410 
SIM04420 
SIM04430 
SIM04440 
SIM04450 
SIM04460 
SIM04470 
SIM04480 
SIM04490 
SIM04500 
SIMa4510 
SIM04520 
SIM04530 
SIM04540 
SIM04550 
SIM04560 
SIM04570 
SIM04580 
SIM04590 
SIM04600 
SIM04610 
SIM04620 
SIM04630 
SIM04640 
SIM04650 
SIM04660 
SIM04670 
SIM04680 
SIM04690 
5IM04700 
SIM04710 
SIM04720 
SIM04730 
SIM04740 
SIM04750 
SIM04760 
SIM04770 
SIM04780 
SIM04790 
SIM04800 
SIM04810 
SIM04820 
SIM04830 
SIM04840 
SIM04850 
SIM04860 
SIM04870 
SIM04880 
SIM04890‘ 
SIM04900 
SIM04910 
SIM04920 
SIM04930 
SIM04940 
SIM04950 
SIM04960 
SIM04970 
SIM04980 
SIM04990 
SIM05000 
SIM05010 
SIM05020 
T0SIM05030 
SIM05040 
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RI  =  (RNEM-RMIN)/DR+. 999999  SIM05050 

IR  ::  RI  SIM05060 

TI  =  (TMAX-T)/DDALFA+. 999999  SIM05070 

IT  =  TI  SIM05080 

IF  (IT.LE.O)GO  TO  314  SIM05090 

IF  (IT.GT.NAD)GO  TO  315  SIM05100 

IF(IR.GT.NRD)GO  TO  316  SIM05110 

C  COUNT  S  DIRECTION,  SAMPLE . 18  SIM05120 

FSN2(1,IT,KR)=FSN2(1,IT,KR)+1.  SIM05130 

GO  TO  319  SIM05140 

C  COUNT  N  DIRECTION,  SAMPLE . 18  SIM05150 

314  IF(IR.LE.0)IR=1  S1M05160 

IF(IR.GT.NRD)IR=NRD  SIM05170 

FWN2(1,IR)=FNN2(1,IR)+1.  S1M05180 

GO  TO  319  SIM05190 

C  COUNT  E  DIRECTION,  SAMPLE . 18  SIM05200 

315  IF(IR.LE.0)IR=1  SIM05210 

IF(IR.GT.NRD)IR=NRD  SIM05220 

FEN2(1,IR)=FEN2C1,IR)+1.  SIM05230 

GO  TO  319  SIM05240 

C  COUNT  N  DIRECTION,  SAMPLE, . 18  SIM05250 

316  FNN2(1,IT,KR)=FNN2(1,IT,KR)+1.  SIM05260 

C  SIM05270 

319  RNEW  =  -99.  SIM05280 

C . 19  SIM05290 

313  P2(4,I2)=  RNEM  SIM05300 

P2(5,I2)=  T  SIM05310 

320  CONTINUE  SIM05320 

C  SIM05330 


C  NON  NEN  MOLECULES  HILL  BE  INTRODUCED 

C  TOTAL  JET  FLUX  INTO  THE  REGION  WAS  DETERMINED  BY  FKI)  MOL. /SEC. 

C  OR  FWP1(I,KR),FWP2(I,KR)  FOR  EACH  SPECIES 

C  NUMBER  OF  MOLECULES  TO  BE  ACTIVATED  PER  DTM  IS  FCDXDTM  PER  CELL. 

C 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

DO  390  I:1,NRD 
C  INPUT  'W*  MOLECULES 

IF(KS.E0.1)G0  TO  322 
ANEW! =F0E1 (1,1, KR, KS-1 ) 

ANEW2=F0E2(1,I,KR,KS-1) 

GO  TO  323 

322  ANEW1=FWP1C1,I,KR)XDTM 
AHEW2=FWP2(1,I,KR)xdTM 

323  NEHl^ANENl 
NEH2-ANEW2 
REM1=ANEW1-NEW1 
CALL  RANDU(P) 

IFCREMl .GT.P)NEW1=NEW1+1 

P.EM2=ANEW2-NEW2 

CALL  RANDU(P) 

IF(REM2.GT.P)NEW2=NEW2+1 
C  ACTIVATE  NEW  INPUT  MOLECULES 

C  TIME, LOCATION  AND  VELOCITY  COMP.  OF  NEW  MOLS.  ARE  RANDOM  FUNCTIONS 

C  SPECIES  1 . 

IFCNEWl.LT.DGO  TO  341 
C 

DO  340  11=1, NEWl 
CALL  RANDU(P) 

atime=pxdtm 

CALL  RANDU(P) 

RSTART=(FLOAT(I-11+P)xDR+RMIN 

VELX=FWP1(2,I,KR) 

VELY=FWP1(3,I,KR) 

TEM=FWP1(4,I,KR) 

C  THERMAL  VELOCITY 

VTER1=SQRT(2.XB0LTZXTEM/SPEC(1,1)) 

C  TEMPERATURE  IS  TRANSFERRED  THROUGH  FWPK 4, I , KR) 
XX=RSTARTxCOS(ALFA)+VELXxATIME 
YY=RSTARTXSIN(ALFA)+VELYXATIME 
RX=SQRT(XXXX2+YYXX2) 

T=ATAN(YY/XX) 


SIM05350 

SIM05360 

SIM05370 

SIM05380 

SIM05390 

SIM05400 

SIM05410 

SIM05420 

SIM05430 

SIM05440 

SIM05450 

SIM05460 

SIM05470 

SIM05480 

SIM05490 

SIM05500 

SIM05510 

SIM05520 

SIM05530 

SIM05540 

SIM05550 

SIM05560 

SIM05570 

SIM05580 

SIM05590 

SIM05600 

SIM05610 

SIM05620 

SIM05630 

SIM05640 

SIM05650 

SIM05660 

SIM05670 

SIM05680 

SIM05690 

SIM05700 

SIM05710 

SIM05720 

SIM05730 

SIM05740 

SIM05750 

SIM05760 
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TWLW'.^  WLWLWl  W 
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IF(RX.LT.RNIN.OR.RX.GT.RNAX)GO  TO  S40 
IF(T.LT.TMIN.OR.T.GT.THAX)GO  TO  340 
C  FOR  MORE  ACCURATE  CALCULATION  SET  THESE  MOLECULES  IN  OUTPUT  FLOHS 
C 

C  DEFINE  THE  VELOCITY  COMPONENTS 
CALL  RANDU(P) 

B*VTER1XSQRT(-AL0G(P)) 

CALL  RANDU(P) 

BBs2.»PI»P 

VELXsVELX-i>BXCOS(BB) 

VELY*VELY+B»SIM<BB) 

CALL  RANDU(P) 

VELZ^^VTERIXSQRTC  -ALOG(P)  ) 

C  DEFINE  THE  NEH  MOLECULE  TO  BE  ACTIVATED 
DO  325  IACT>1,NM0L1 
IF(P1(4,IACT).EQ.-99.)G0  TO  326 

325  CONTINUE 

IF  THERE  IS  NO  ROOM  FOR  ADDITIONAL  MOLECULE  PRINT  'ALARM* 
IFdACT.GE.NMOLDGO  TO  3004 

326  P1(4,IACT)-RX 
P1(5,IACT)=T 
P1(1,IACT)=VELX 
P1C2,IACT)=VELY 
P1(3,IACT)=VELZ 

340  CONTINUE 

' REPEAT ' PROCEDURE ‘ FOR ‘ SPEcIeS  *2 . 

341  IF(NEW2.LT.1)G0  TO  351 
DO  350  12-1, NEM2 
CALL  RANDU(P) 

ATIME=P»DTM 
CALL  RANOU(P) 

RSTART=RMIN+(FL0AT(I-1)+P)»DR 
VELX=FWP2(2,I,KR) 

VELY=FWP2(3,I,KR) 

TEM=FWP2(4,I,KR) 

THERMAL  VELOCITY 

VTER2»SQRT(2.»B0LT2*TEM/SPEC(2,1)) 
XX=RSTARTXCOSCALFA)+VELXXATIME 
YY=RSTARTXSIN(ALFA)+VELYXATIME 
RX=S0RT(XXX»2+YYX»2) 

T=ATAN(YY/XX) 

IF(RX.LT.RMIN.OR.RX.GT.RMAX)GO  TO  350 
IF(T.LT.TMIN.OR.T.GT.TMAX)GO  TO  350 
FOR  MORE  ACCURATE  RESULTS  SET  THESE  MOLECULES  IN  OUTPUT  FLOWS 

DEFINE  VELOCITY  COMPONENTS 
CALL  RANDU(P) 

B»VTER2XS(}RT(-AL0G(P)) 

CALL  RANDU(P) 

BB=2.XP1XP 
VELX=VELX+BXC0S(BB) 

VELY=VELY+BXSIN(BB) 

CALL  RANDU(P) 

VEL2=VTER2XSQRT(-AL0G<P) ) 

FIND  A  NEW  MOLECULE  TO  BE  ACTIVATED 
DO  345  IACT=1,NM0L2 
IF  (P2(4,IACT).EQ.-99.)G0  TO  346 

345  CONTINUE 

IF  THERE  IS  NO  PLACE  THEM  PRINT  ALARM 
IF(IACT.EQ.NM0L2)G0  TO  3004 

346  P2(4,IACT)=RX 
P2(5,IACT)=T 
P2(1,IACT)=VELX 
P2(2,IACT)=VELY 
P2(3,IACT)=VELZ 

350  CONTINUE 

REPEAT  PROCEDURE  FOR  SPEC-3 

INPUT  E  MOLECULES . 


SIM05770 

SIM05780 

SIM0S79O 

SIM05800 

SIM05810 

SIN05820 

SIM0S830 

SIM0S840 

SIM05850 

SIM05860 

SIM05870 

SIM05880 

SIM0S890 

SIM05900 

SIM05910 

SIM05920 

SIM05930 

SIM0S940 

SIM05950 

SIM05960 

SIM05970 

SIM05980 

SIM05990 

SIM06000 

SIM06010 

SIM06020 

SIM06030 

SIM06040 

SIM06050 

SIM06060 

SIM06070 

SIM06080 

SIM06090 

SIM06100 

SIM06110 

SIM06120 

SIM06130 

SIM06140 

SIM06150 

SIM06160 

SIM06170 

SIM06180 

SIM06190 

SIM06200 

SIM06210 

SIM06220 

SIM06230 

SIM06240 

SIM06250 

SIM06260 

SIM06Z70 

SIM06280 

SIM06290 

SIM06300 

SIM06310 

SIN06320 

SIM06330 

SIM06340 

SIM06350 

SIM06360 

SIM06370 

SIM06380 

SIM06390 

SIM06400 

SIM06410 

SIM06420 

SIM06430 

SIM06440 

SIM06450 

SIM06460 

.SIM06470 

SIM06480 
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351  IF  (TMIN.LE.O. )G0  TO  390 


ANEHlsFOHl ( 1 , 1 , KR»  KS^l ) 


NEHl-ANEWl 
REN-ANEHl-NEHl 
CALL  RANOU(P) 
IF(P.LT.REN)NEN1*NEH1+1 

ANEN2-F0N2( 1 , I . KR, KS^l ) 
NEM2>ANEN2 
REM=ANEN2-NEM2 
CALL  RANDU(P) 

1F(  P .  LT .  REM)NEH2-NEH24^1 


IF(NEN1.LE.O)60  TO  370 

DO  362  I1’‘1,NEM1 
CALL  RAN0U(P> 

T=TMIN+DDALFA*P 
CALL  RANOU(P) 

RNEMs  RMI  N-f  (  FL  OAT  ( I  - 1 ) -t-P  )  XDR 
CALL  RANDU(P) 

TEM=F0N1 C  4 , 1 , KR. KS+1 ) 

VTERl =SQRT( 2 . xBOLTZXTEM/SPECt 1 , 1 ) ) 

B=VTER1XS(}RT(-AL0G(P)  ) 

CALL  RANOU(P) 

BB=2.XPIXP 

VELX^FOHl  (2, 1 ,  KR,  KS-t^l  )>BXC0S(BB) 

VELY-FOMl ( 3, 1 . KR, KS-l-1  )-i-BXSIN( BB) 

CALL  RANDU(P) 

VELZ=VTER1»SQRT(-ALQG(P)) 

FIND  A  NEM  MOLECULE  TO  BE  ACTIVATED 
DO  364  lACT^l.NMOLl 
IF(P1(4,IACT).EQ.-99.)G0  TO  366 
364  CONTINUE 

IF  THERE  IS  NO  ROOM  FOR  ADDITIONA*  MOLECULE  THEN  PRINT  'ALARM* 
IFdACT.EQ.NMOLDGO  TO  3004 
366  P1(4,IACT)-RNEH 
P1<5,IACT)=T 
P1<1,IACT)=VELX 
P1(2,IACT)*VELY 
362  P1C3,IACT)=VELZ 
370  CONTINUE 

REPEAT  FOR  SPECIES  2 

IF(NEN2.LE.O)GO  TO  390 
DO  382  12=1, NEM2 
CALL  RANDU(P) 

T=TMIN+DDALFAKP 
CALL  RANDU(P) 

RNEW=RMIN+( FLOAT! I-1)+P)KDR 


S1N06490 

SIM06500 

SIM06510 

SIN06520 

SIM06530 

SIM06540 

SIM06550 

SIM06560 

SIM06S70 

SIM06580 

SIM06590 

SIM06600 

SIM06610 

SIM06620 

SIM06630 

SIM06640 

SIM06650 

SXM06660 

SIM06670 

SIM06680 

SIM06690 

SIM06700 

SIM06710 

SIM06720 

SIM06730 

SIM06740 

SIM06750 

SIM06760 

SIM06770 

SIM06780 

SIM06790 

SIM06800 

SIM06810 

SIM06820 

SIM06830 

SIM06840 

SIM06850 

SIM06860 

SIM06870 

SIM06880 

SIM06890 

SIM06900 

SIM06910 

SIM06920 

SIM06930 

SIM06940 

SIM069S0 

SIM06960 

SIM06970 

SIM06980 

SIM06990 

SIM07000 

SIM070I0 

SIM07020 

SIM07030 


CALL  RANDU(P) 

TEMsF0W2{4,I,KR,KS+l) 

VTER2=SQRT(  2 . xBOLTZxTEM/SPECCZ, 1 ) ) 
B=VTER2*S0RT(-AL0G(P)) 

CALL  RANDU(P) 

BB=2.»PI*P 

VELX=F0W2C2,I,KR,KS+1)+BkC0S(BB) 
VELY=F0W2( 3, I , KR, KS+1 )+BKCOSCBB) 
CALL  RANDU(P) 
VELZ=VTER2»S0RT(-AL0G(P) ) 

C  FIND  A  NEH  MOLECULE  TO  BE  ACTIVATED 
DO  384  IACT=1,NM0L2 


SIM07040 

SIM07050 

SIM07060 

SIM07070 

SIM07080 

SIM07090 

SIM07100 

SIM07110 

SIM07120 

SIM07130 

SIM07140 

SIM07150 


IF(P2(4,IACT).EQ.-99)G0  TO  386 
384  CONTINUE 

C  IF  THERE  IS  NO  ROOM  FOR  ADDITIONAL  MOLECULE  THEN  PRINT  ALARM 
IF(IACT.EQ.NM0L2)  GO  TO  3004 
386  P2(4,IACT)=RNEM 


SIM07160 

SIM07170 

SIM07180 

SIM07190 

SIM07200 


'  ■.* 
•-.  V  O  •_ 


OOO 


P2<5,IACT)»T 
P2(1.IACT)’:VELX 
P2(2,IACT)*VELY 
382  P2(3.IACT)-VELZ 


REPEAT  FOR  SPEC-3 
390  CONTINUE 


INPUT  'S'  MOLECULES 
DO  460  J-1,NAD 
IF(KR.LT.2)60  TO  460 
ANEM1-FNN1(1,J.KR-1) 

IFCANEHl.LE. .00001)60  TO  420 
NEN1»ANEM1 
REM-ANEWl-NEMl 
CALL  RANDU(P) 
IF(P.LT.REN)NEMlsNEWl4l 
DO  402  I1»1,NEM1 
CALL  RANDU(P) 

T=TMIN+DDALFA»( P+FL0AT( J-1 ) ) 

CALL  RANDU(P) 

RNENsRMIN+PXDR 
CALL  RANOU(P) 

TEM*FNN1(4, J,KR-1) 
VTER1=SQRT(2.»B0LTZXTEM/SPEC(1,1)) 
B=VTER1»SQRT(-AL0G(P) ) 

CALL  RANDU(P) 

BB*2 • xPixp 

VELX=FNM1(2.J.KR-1)+B»C0S(BB) 
VELY-FNN1(3.  J,KR-l)-fBXSIN(BB) 

CALL  RANDU(P) 
VEL2=VTERU«SQRT(-AL0G<P)) 

C  FIND  A  NEM  MOLECULE  TO  BE  ACTIVATED 
DO  404  IACT-1,NM0L1 
IF(P1(4,IACT) .EQ.-99. )G0  TO  406 
404  CONTINUE 

IFCIACT.EQ.NMOLDGO  TO  3004 
406  P1(4,IACT)=RNEH 
P1(5,IACT)»T 
Pl(l,IACT)sVELX 
P1(2,IACT)»VELY 
402  P1C3,IACT)=VELZ 
420  CONTINUE 
C  REPEAT  FOR  SPECIES  2 

ANEM2=FNN2(1, J,KR-1) 

IF(ANEM2.LT.  .OOOODGO  TO  440 

NEN2=ANEH2 

REM=ANEW2-NEW2 

CALL  RANOU(P) 

IF(P.LT.REM)NEW2*NEW2+1 

DO  422  12=1, NEM2 

CALL  RANDU(P) 

T=TMIN+DDALFAKCP+FL0AT(J-1)) 

CALL  RANDU(P) 

RNEM=RMIN->-PxDR 
CALL  RANDU(P) 

TEM=FNN2(4, J,KR-1) 

VTER2=S0RT<2.*B0LTZ*TEM/SPEC(2,1)) 

B=VTER2XS0RT<-AL0G(P)) 

CALL  RANDU(P) 

BB=2.»PIKP 

VELX=FNN2(2,J,KR-1)>BXC0S(BB) 
VELY=FNN2(3, J,KR-1)+BKSIN(BB) 

CALL  RANDU(P) 
VELZ=VTER2KSQRT(-AL0G(P)) 

C  FIND  A  NEM  MOLECULE  TO  BE  ACTIVATED 
DO  424  IACT=1,NM0L2 
IF(P2(4,IACT).EQ. -99)00  TO  426 
424  CONTINUE 

IF(IACT.EQ.NM0L2)G0  TO  3004 
426  P2C4,IACT)=RNEM 
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P2(5,IACT)»T 
P2(1,1ACT)=VELX 
P2(2,IACT)*VELY 
422  P2(3>ZACT}sVELZ 
440  CONTINUE 


REPEAT  FOR  SPEC-3 


INPUT  'N*  MOLECULES . 

SPEC-1 

ANEM1-FSN1(1,J.KR+1) 

IF(ANEM1.LE. .00001)60  TO  450 

NEH1>ANEM1 

REN>ANEH1-NEW1 

CALL  RANDU(P) 

IF  (P.LT.REM)NEN1>NEM1+1 
DO  442  I1«1,NEM1 
CALL  RANDU(P) 

T=TMIN  ♦DDALFA*(P+FL0AT(J-1)) 

CALL  RANDU(P) 

RNEM*RMAX-PXDR 
CALL  RANDU(P) 

TEM=FSN1(4,J,KR+1) 

VTER1=SQRT(2.»B0LTZ«TEM/SPEC(1,1)) 

B=VTER1*SQRT(-AL06(P)) 

CALL  RANOU(P) 

BB-2  XPIXP 

VELXsFSNl(2.J.KR-*-l)-«-BXCQS(BB) 

VELY=FSNl(3,J,KR-fl}-*-Bi(SIN(BB) 

CALL  RANDU(P) 
VELZ=VTER1*S0RT<-AL0G(P)) 

FIND  A  NEM  MOLECULE  TO  BE  ACTIVATED 
DO  444  IACT-1,NM0L1 
IF(Pl(4,IACT).EQ.-99.)60  TO  446 
444  CONTINUE 

IFCIACT.EQ.NMOLDGO  TO  3004 
446  P1(4,IACT)=RNEW 
Pl<5,IACT)sT 
P1(1,IACT)=VELX 
P1C2,IACT)»VELY 
442  P1C3,IACT)=VELZ 
450  CONTINUE 


INPUT  N  MOLECULES  SPEC-2 
ANEW2»FSN2(1,J,KR+1) 

IF(ANEH2.LE. .00001)60  TO  460 
NEH2-ANEM2 
REM=ANEW2-NEM2 
CALL  RANDU(P) 

IFC  P . LT . REM)NEM2=NEW2+1 
DO  452  12-1, NEW2 
CALL  RANDU(P) 

T=TMIN+DDALFA»<P+FL0AT<J-1)) 

CALL  RANDU(P) 

RNEM-RMAX-PXDR 
CALL  RANDU(P) 

TEM=FSN2(4, J,KR+1) 

VTER2  =  S0RT<  2 . »BOLT2kTEM/SPEC( 2, 1 ) ) 
B=VTER2»S0RT<-AL0G(P)) 

CALL  RANDU(P) 

BB=2»PI*P 

VELX=FSN2(2, J,KR+l)+BXCOS<BB) 
VELYsFSN2(3, J,KR+1)+B»3IN(BB) 

CALL  RANDU(P) 
VELZaVTER2*SQRT(-AL0G<P)) 

C  FIND  A  NEM  MOLECULE  TO  BE  ACTIVATED 
DO  454  1ACT=1,NM0L2 
IF(P2(4,IACT).EQ. -99)60  TO  456 
454  CONTINUE 

IF(IACT.EQ.NM0L2)G0  TO  3004 
456  P2(4,IACT)-RNEH 
P2(5,1ACT)=T 
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P2C1.IACT)*VELX 
P2(2,IACT)»VELY 
452  P2(3,IACT)sVELZ 
460  CONTINUE 
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REARRANGE  MOLECULES  IN  THEIR  CELLS 
INITIALIZATION 
5000  CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
KIPsNNOLl-»^NMOL2-«-NMOL3 
00  1001  KAOO^l.KIP 
1001  IP(KADD)-0 
KADD-0 
N1A=0 
N2AsO 
N3A»0 

N1A,N2A,N3A  ARE  THE  NUMBER  OF  ACTIVE  MOLECULES  (COUNTED  NEXT) 

DO  1100  I»1,NRD 
DO  1100  J=l,NAO 

DO  1005  K»l,4 
1005  IC(K,I.J)-0 

SET  SPECIES  1 
N1C«0 
N2C=0 
N3CsO 
NTC=0 

DO  1020  K1-1,NM0L1 
IF(P1(4,K1).EQ.-99.)G0  TO  1020 
RL0C-ABS(C(19,I, J)-P1(4,K1)) 

TLOC=ABS(C(20.I. J)-P1(5,K1)) 

IF(RL0C.GT.DRX.5)G0  TO  1020 
1F(TL0C.GT.DDALFAX.5)G0  TO  1020 

N1C=N1C+1 

NTC=NTC+1 

N1A=N1A+1 

IFCNTC.EQ.l)IC(4,I,J)sKADD 

IC(1,I,J)=IC<1,I,J)+1 

KADD3KADD41 

IP(KADD)=K1 

1020  CONTINUE 

SET  SPEC-2  MOLECULES . 

DO  1040  K2=1,NM0L2 
IF(P2(4,K2) .EQ.-99. )GO  TO  1040 
RL0C=ABS<C(19,I,J)-P2(4,K2)) 

TLOC=ABS(C(20,I,J)-P2(5,K2)) 

IF(RL0C.GT.DRX.5)G0  TO  1040 
IF(TL0C.GT.DDALFAX.5)G0  TO  1040 

N2A:N2A-*-l 

NTC=NTC-M 

IF(NTC.EQ.niC(4,I,J)=KADD 
ICC2,I,J)  =  IC(2,I,  J)  +  l 
KADD^KADD-^l 
IP(KADD)=K2 
1040  CONTINUE 
1100  CONTINUE 

MRITE(6,1021)N1A,N2A 

1021  FORMATC  NUMBER  OF  ACTIVE  MOLECULES  SPEC1=  M5,  •  SPEC2=M5/) 
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DO  999  I-1»NR0 
00  999  J«1»NAD 
999  NC0L(Z.J)»0 
C 

DO  900  I>1,NR0 
DO  900  J*1,NA0 
C 

NM(l)slC(l,I,J) 

NM(2)-IC(2.I»J) 

NH(3)-IC(3»I,J) 

NTOTsNMC  1  )4>NM(2)+NM(  3) 

DO  900  L»l»3 
DO  900  M>1,3 

C  NUMBER  OF  SPECIES  IN  THE  PROGRAM  IS  2  (3) 

KVs(L-l)*3+M 

KT*KV+9 

C  KV,KT  ARE  THE  ADDRESSES  OF  RELATIVE  VELOCITY  AND  COLLISION  TIMERS 
920  IF(C(KT,I.J).GE.TIME)  GO  TO  900 
C  C(KT.I,J)IS  THE  INTEGRATED  TIME  FOR  L-M  COLLISION 
KSEL-0 
KREJ-0 

1F(NM(L).GT.1.AND.NM(M).GT.1)G0  TO  912 
C  NO  COLLISIONS  ARE  CALCULATED  IF  THERE  ARE  NO  MOLECULES 

911  CCKT.I,J)=C(KT,I.J)+DTM 
GO  TO  900 

C  SELECT  NOM  THE  MOLECULES  FOR  COLLISION 

912  IF  (KSEL.GE.IOO)GO  TO  911 
CALL  RANOU(P) 

MOLlsPXNM<L)+. 999999 
lF(MOLl.Eq.O)MOLl-l 

C 

CALL  RANDU(P) 

M0L2*P*NM(M)+. 999999 
IF(MOL2.EQ.O)MOL2-1 
C 

KSEL^KSEL+l 

CHECK  IF  THE  SAME  MOLECULE  HAS  BEEN  SELECTED  TWICE 
IF(L.EQ.N.AN0.M0L1.EQ.M0L2)60  TO  912 


FIND  THE  ACTUAL  ADDRESSES  OF  THE  SELECTED  MOLECULES 
IFCL.EQ.DKlsO 
IF(L.EQ.2)K1sNM(1) 

IF(L.Eq.3)Kl3NM(l)-*-NM(2) 

IF(M.EQ.1)K2=0 
lF(M.Eq.2)K2=NM(l) 

IF(M.EQ.3)K2-NM(1)4-NM(2) 

KA01sM0Ll-fKl-MC(9,I,  J) 

KAD2=M0L2+K24^IC(4,I,  J) 

KAD1,KAD2  ARE  THE  LOCATION  OF  SELECTED  MOLECULES  IN  IP(  ) 
MADl-IP(KADl) 

MAD2=1P(KA02) 

MAD1,MAD2  ARE  THE  ACTUAL  ADDRESSES  OF  THE  SELECTED  MOLECULES  (THE 
INDICATION  OF  WHAT  SPECIES  THEY  ARE  HAS  BEEN  DEFINED  HERE) 

DO  930  N=l,3 

IF(L.EQ.l)VNlsPl(N,MADl) 

IF(L.EQ.2)VN1=P2(N,MAD1) 
lF(L.EQ.3)VNlsP3(N,MA01} 

IF(M,E0.1)VN2=P1(N,MAD2) 

IF(M.E0.2)VN2»P2CN,MAD2) 

IF(M.EQ.3)VN2=P3(N,MAD2) 

930  VRC(N)*VN1-VN2 
C  VRC(3)  COrlTAIN  THE  THREE  RELATIVE  VELOCITY  COMPONENTS 
VR=S0RT(VRC<1)»VRC(1)+VRCC2)»VRC(2)+VRC(3)»VRC(3)) 

C  VR  IS  THE  RELATIVE  SPEED  IN  A  SPECIFIC  COLLISION 
IF(C(KV,I, J) .LT.VR)C<KV,I,J)=VR 

C  LAST  STATEMENT  RESETS  THE  MAXIMUM  RELATIVE  VELOCITY  FOR  FURTHER 
C  CALCULATIONS 
C 

IF(KREJ.6T.100)GO  TO  911 


SIM09S70 

SIM09380 

SIM09390 

SIM09400 

SIM09410 

SIM09420 

SIM09430 

SIM09440 

SIM09450 

SIM09460 

SIM09470 

SIM09480 

SIM09490 

SIM09500 

'SIM09510 

SIM09520 

SIM095S0 

SIM09540 

SIM09550 

SIM09560 

SIM09570 

SIM09580 

SIM09590 

SIM09600 

SIM09610 

SIM09620 

SIM09630 

SIM09640 

SIM09650 

SIM09660 

SIM09670 

S1M09680 

SIM09690 

S1M09700 

SIM09710 

SIM09720 

SIM09730 

SIM09740 

SIM09750 

.SIM09760 

SIM09770 

SIM09780 

SIM09790 

SIM09800 

SIM09810 

SIM09820 

SIM09830 

S1M09840 

SIM09850 

SIM09860 

SIM09870 

SIM09880 

SIM09890 

SIM09900 

SIM09910 

SIM09920 

SIM09930 

SIM09940 

SIM09950 

SIM09960 

SIM09970 

SIM09980 

SIM09990 

SIMIOOOO 

SIMlOOlO 

SIM10020 

SIM10030 

SIM10040 

SIM10050 

SIM10060 

SIM10070 

SIM10080 


«'>o  o  oo  no  no  oo  n  o  f>  n«)  no  oooooo 


C 


CALL  RANDU(P) 

AVR>VR/C(KV,I,J) 

KREJ-KREJ4-1 
IF(AVR.LT.P)GO  TO  912 

LAST  STATEMENT  REJECTS  THE  CALCULATED  COLLISION 
NOH  A  SPECIES  L-H  COLLISION  HAS  BEEN  SELECTED 

CALCULATE  NOW  THE  PROBABILITY  THAT  SUCH  A  COLLISION  MILL  BE  COUNTED 
FOR  THE  L  AND  N  SPECIES  RESPECTIVELY 
LP=1 
LM-1 

CALL  RANDU(P) 

ANM3FL0AT( NM( L ) )/FLOAT( NN(M) ) 

IF(ANM.6T.1.)G0  TO  950 
IF(ANM.LT.P)NP:iO 
GO  TO  955 


950  ANMsl./ANM 

IF(ANM.LT.P)LP»0 


955  CXS=PIX(SPEC(L,2)+SPEC(M,2))*»2/«. 

ALP=LP 

ALM=LM 

V0LUME-REG(9,KR,KS)XV0L(I) 

. . .  ••••  . . .  CHECK 

DNL sFLOAT ( NM( L )  ^VOLUME . 

DNM=FL0AT(NM(M))/V0LUME 
USE  EQ.10.3 

TOC(L,M)=ALP/(CXSXDNM»VR*NM(L))+ALM/(CXS*DNL»VR»NM(M)) 

SET  THIS  VALUE  INTO  C(KT,I,J) 

C(KT,I,J)=C<KT,I,J)+TOC<L,M) 


SAMPLE  THIS  COLLISION 

NC0L(I«J)-NC0L(I,J)->-l 

FIND  RELATIVE  MASSES 
CALL  RANDU(P) 
BB=1.-2.XP 
AA=S0RT(1.-BBXBB) 
VRC(1)-BBXVR 
CALL  RANOU(P) 
BB=2.XPIXP 

VRC(2)=AAXC0S(BB)XVR 

VRC(3):AAXSIN(BB)XVR 

FIND  RELATIVE  MASSES 

SM=SPEC(L,1)+SPEC(M,1) 

RML=SPEC(L,1)/SM 

RMM=SPECCM,1)/SM 
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CALCULATE  HERE  THE  ACTUAL  ADDRESSES  OF  COLLIDING  MOLECULES  AND  SET 
THEIR  NEH  VELOCITIES  (NEM  VELOCITY  COMPONENTS  ARE  ADDED  TO  THE 
VELOCITY  OF  THE  CENTER  OF  MASS  VCCM) 

DO  960  N=l,3 
IF(L.E0.1)V1=P1(M,MAD1) 

IF(L,EQ.2)V1=P2<N,MAD1) 

IFCL.EQ.S)V1=P3(N,MAD1) 

1F(M.EQ.1)V2=P1(N,MAD2) 

IF(M.E0.2)V2=P2(N,MAD2) 

IF(M,EQ.3)V2=P3(N,MAD2) 

VCCM=RMLXV1+RMMXV2 
VCCMl =VCCM+VRC( N ) XRMM 
VCCM2=VCCM-VRC( N ) XRML 
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C  CHANGE  IN  VELOCITY  IS  INPUTEO  ONLY  IF  PROBABILITY  OF  COLLISION. 6T.1  SIN10810 

IFCLP.NE.DGO  TO  961  SIN10820 

IF(L.EQ.l)Pl(N,MADl)sVCCMl  SIM10830 

IF(L.EQ.2)P2(N,MAD1)*VCCM1  SIM10840 

C  IF(L.EQ.3)P3(N,MAD1)*VCCM1  SIM10850 

C  SIM10860 

961  IFCLN.NE.DGO  TO  960  SIM10870 

IF(N.EQ.1)P1(N«MA02)«VCCM2  S1M10880 

IF(M.Eq.2)P2(N,HAD2)>VCCM2  SIM10890 

C  lF(H.EQ.3)P3(N,MAD2)sVCCM2  SIM10900 

C  SIM10910 

960  CONTINUE  SIM10920 

GO  TO  920  SIM10930 

900  CONTINUE  SIM10940 

MRITE(6,4545)T0C(1.1).T0C(1,2)»T0C(2,1),T0C(2,2)»TIME  SIM10950 

4545  FORHATC  TINE* .5E10.3)  SIN10960 

C  SIM10970 

CxxkxkkxxxxxxxxxxiixeND  OF  COLLISIONSxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxSIM10980 


C  NON  NEM  temperatures  MAY  BE  CALCULATED  IN  EACH  CELL 
C  AVERAGE  VELOCITY  IN  CELLS 
C 

MRITE(6,4547) 

4547  FORMATCO  I  J  NIE  N2E  TEMPRl  TEMPR2 

1  VAVYl  VAVX2  VAVY2  NCOL',/) 

C 

DO  1110  I-l.NRO 
DO  1110  J31,NA0 
NlE^ICd.I,  J) 

N2E-IC(2,I,J) 

KAD1»1C(4,I> J) 

KAD2=IC(4,1.J)-^N1E 
IFCNIE.LT.DGO  TO  1112 
C 

VX1=0. 

VX2=0. 

VY1»0. 

VY2=0. 

C 

DO  1111  IM-1,N1E 

MADl^KADUIM 

lAD-IP(MADl) 

VX1=VX1+P1(1,1AD) 

1111  VY1=VY1+P1(2,IAD) 

1112  IF(N2E.LE.O)GO  TO  1110 
C 

DO  1115  IMsl,N2E 

MAD2=MAD1>IM 

IAD-IP(MAD2) 

VX2=VX2+P2C1,IAD) 

1115  VY2=VY2*P2C2,IAD) 

C 

C  THE  AVERAGE  VELOCITIES  IN  THE  CELL  HILL  RESULT- 
VAVX1=VX1/N1E 
VAVY1»VY1/N1E 
VAVX2=VX2/N2E 
VAVY2=VY2/N2E 
IFCI.NE.DGO  TO  1116 
C 

FSN1(2,J,KR)=VAVX1+FSN1(2,J,KR) 

FSH1(3, J,KR)sVAVYl+FSNl(3, J,KR) 

FSH2(2, J.KR)=VAVX2+FSN2(2,J,KR) 
FSN2(3,J.KR)-VAVY2-i-FSN2(3,J,KR) 

C 

1116  IF(I.NE.NRD)GO  TO  1117 
FNN1(2,J,KR)=VAVX1+FNN1(2,J,KR) 

FNN1C3, J,KR)*VAVY1+FNN1C3,J,KR) 

FNN2( 2 , J . KR ) =VAVX2*FNN2( 2,J,KR) 

FNN2 ( 3 , J , KR ) =VAVY2+FMN2<  3 , J , KR ) 

C 

1117  IFCJ.NE.DGO  TO  1118 
FEN1(2,I)»VAVX1+FEN1C2,I) 

FEN1<3,1)=VAVY1+FEN1(3,I) 


SIM10990 
SIMllOOO 
SIMllOlO 
SIM11020 
VAVXl  SIM11030 
SIM11040 
SIM11050 
SIM11060 
SIM11070 
SIM11080 
SIM11090 
SIMlllOO 
SIMllllO 
SIN11120 
SIM11130 
SIM11140 
SIM11150 
SIM11160 
SIM11170 
SIM11180 
SIM11190 
SIM11200 
SIM11210 
SIM11220 
SIM11230 
SIM11240 
S1M11250 
S1M11260 
SIM11270 
SIM11280 
SIM11290 
SIM11300 
S1M11310 
SIM11320 
SIM11330 
SIM11340 
SIM11350 
SIM11360 
SIM11370 
SIM11380 
SIM11390 
SIM11400 
SIM11410 
SIM11420 
SIM11430 
SIM11440 
SIM11450 
SIM11460 
S1M11470 
SIM11480 
SIM11490 
SIM11500 
SIM11510 
S1M11520 
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FEN2( 2, 1 ) >VAVX2+FEN2 (2*1) 

FEN2( 3, 1 >  *VAVY2>FEN2( 3. 1 ) 

C 

1118  IF(J.NE.NAD)60  TO  1119 
FHN1(2>I)=VAVX1«FNN1(2,I) 
FMNl(3,l)-VAVYl-i-FNNl(3,I) 

FMN2(  2 , 1 )  =VAVX2+FMtl2(  2  >  I ) 

FNN  2  (  3 . 1 )  >  VA  VY2-^  FNN2  (3.1) 

1119  CONTINUE 

THERNAL  VELOCITIES  AND  TEMPERATURES 
ENR61-0. 

ENR62-0. 

IF(N1E.LT.1)60  T01131 
DO  1130  IM^l.NlE 
MAD1»KAD1-»^IN 
IA0>IP(MA01) 

CXlsPl(l,IAD)-VAVXl 
CY1=P1(2,IAD)-VAVY1 
CZ1-P1(3.IAD) 

1130  ENRGl  - ENRGl -t^CXl XCXl -t-CYl XCYl-f CZl xCZl 

1131  IF(N2E.LT.1)GQ  TO  1150 
DO  1140  IM=1,N2E 
MAD2=KAD2-MN 
lAD^IPCMAOZ) 

CX2=P2(1,IA0)-VAVX2 
CY2=P2(2,IAD)-VAVY2 
CZ2sP2(3>lA0) 


1140  ENRG2S ENRG2+CX2XCX2+CY2XCY2+CZ2XCZ2 
C 

1150  TEMPR1=(ENRG1/M1E)XSPEC(1,1)/(3.»B0LTZ) 
TEMPR2=(ENRG2/N2E)»SPEC(2,1)/(3.*B0LTZ) 
C 

IFd.NE.DGO  TO  1151 

FSN1(4, J,KR)=TEMPR1+FSN1(4,J.KR) 

FSN2(4,J,KR)=TEMPR2+FSN2(4,J,KR) 

'’1151  IF<I.NE.NRD)GO  TO  1152 

FMNl<4,J,KR)sTEMPRl+FmU(4,J,KR) 

FNN2 ( 4 , J , KR ) =TEMPR2+FNN2 ( 4 . J . KR ) 

C 

1152  IFCJ.NE.DGO  TO  1153 

FWN1(4,I)=TEMPR1+FWN1(4,I) 

FWN2( 4 , I ) =TEMPR2+FHN2( 4 , I ) 

*'1153  IF(J.NE.NAD)GO  TO  1154 


SIM11530 

SIM11540 

SIM11550 

SIM11560 

S1M11570 

SIM11580 

SIM11590 

SIM11600 

S1M11610 

S1M11620 

SIM11630 

SIM11640 

SIM11650 

SIM11660 

SIM11670 

S1M11680 

SIM11690 

S1M11700 

SIM11710 

S1M11720 

SIM11730 

SIM11740 

SIM11750 

S1M11760 

SIM11770 

SIM11780 

SIM11790 

SIM11800 

SIM11810 

SIM11820 

SIM11830 

SIM11840 

SIM11850 

SIM11860 

SIM11870 

SIM11880 

SIM11890 

S1M11900 

SIM11910 

S1M11920 

S1M11930 

SIM11940 

SIM11950 

S1M11960 

S1M11970 

SIM11980 

SIM11990 

simzooo 


FEN1(4,I)=TEMPR1+FEN1(4,I)  SIM12010 

FEN2(4,I)=TEMPR2+FEN2(4,I)  S1M12020 

1154  CONTINUE  SIM12030 

C  SIM12040 

WRITE(6,4548)I.J.N1E.N2E.TEMPR1,TEMPR2.VAVX1.VAVY1.VAVX2.VAVY2.NCOSIM12050 
ILd.J)  SIM12060 

4548  FORMATC*  ' ,4I5,6F12 . 3, 17 )  SIM12070 

C  SIM12Q8Q 


1110  CONTINUE 


SIM12090 


C 

6000  CONTINUE 

C  CALCULATE  AVERAGED  PARAMETERS,  WEIGHTED  BY  DFI 
DO  6010  1=1, NRD 
DO  6010  KPAR=1,4 

F0W1(KPAR,I,KR,KS)=FWN1(KPAR,I)/(NIS) 

F0W2(KPAR.I,KR,KS)=FWN2(KPAR,I)/(NIS) 

F0E1(KPAR,I,KR,KS)=FEN1(KPAR,I)/(NIS) 

6010  F0E2( XPAR, I , KR, KS) =FEN2(KPAR, I )/(NIS) 

C 

DO  6020  J=1,NAD 
DO  6020  KPAR=1,4 

FSNKKPAR,  J,KR)  =  FSN1(KPAR,  J.KR)/NIS 
FSN2(KPAR, J,KR)=FSN2(KPAR, J,KR)/NIS 
FNN1(KPAR,J,KR)=FNN1(KPAR,J,KR)/NIS 


SIM12100 

SIM12110 

SIM12120 

SIM12130 

SIM12140 

SIM12150 

SIM12160 

SIM12170 

SIM12180 

SIM12190 

SIM12200 

SIM12210 

SIM12220 

SIM12230 

SIM12240 


6020  FNN2(KPAR,J,KR)sFNN2(KPAR,J,KR)/NIS  SIM12250 
C  SIM12260 
C  TO  PREPARE  FLOWS  FOR  THE  NEXT  REGION  OR  SECTOR,  DIVIDE  FONO  BY  L0CALSIM12270 


C  DFI  AND  MULTIPLY  BY  NEXT  DFI  WHEN  STARTING  NEW  REGION  SIM12280 
C  SIM12290 
C  CALCULATE  HERE  THE  MEAN  FREE  PATH  AND  STORE  INTO  REG(,,)  SIM12300 
C  SIM12310 
C  STOP  PROGRAM  IF  FLOW  BECOMES  COLLISIONLESS  OR  NUMBER  DENSITY  IS  EQUAL  SIM12S20 
C  TO  THE  AMBIENT  NUMBER  DENSITY.  STORE  F0E1,F0E2  IN  A  SEPARATE  FILE  TO  SIM12330 
C  BE  USED  IN  A  DIFFERENT  PROGRAM.  SIM12340 


C  PRINT  AVERAGED  RESULTS 
WRITE(6,6029) 

6029  FORMATC  ',✓///) 

WRITE(6,6030)NIS.KR 


SIM12340 

SIM12350 

SIM12S60 

SIM12370 

SIM12380 

SIM12390 


6030  FORMATC  AVERAGED  OUTPUT  FLOWS  AFTER*, ISC  TIME  INCREMENTS  IN  REGISIN12400 


ION*, IS) 

WRITE(6,60S1) 

6031  FORMATC *0  I  FEN1(1,I)  FEN2(1,I)  FWN1(1,I)  FWN2 

1) 

DO  6033  I»1.NRD 

WRITE(6.6032)I,FEN1(1,1),FEN2(1,I),FWN1(1,I),FWN2(1,I) 

6032  FORMATC  •,I5,4F13.3) 

6033  CONTINUE 
WRITE  (6,6034) 

6034  FORMATCO  J  FNN1(1,J,KR)  FNN2(1,J,KR)  FSN1(1,J,KR) 
12(1, J,KR)*) 

DO  6036  J^l.NAD 


SIM12410 
SIM12420 
FWN2(1,I)*SIM12430 
Sir>ll2440 
SIM12450 
SIM12460 
SIM12470 
SIM12480 
SIM12490 
,KR)  FSNSIM12500 
SIM125I0 
SIM12520 


WRITE(6,6037)J,FNN1(1,J,KR),FNN2(1,J,KR),FSN1(1,J,KR),FSN2(1,J,KR)SIM12530 


6037  FORMATC  •,I5,4F15.5)  SIM12540 

6036  CONTINUE  SIM12550 

C  SIM12560 

cxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx  SIM12570 
C  SIM12580 

C  START  A  NEW  REGION  SIM12590 

C  IFCKR.EQ.IO)  GO  TO  7000  SIM12600 

C  KR^KR+l  SIM12610 

C  GO  TO  2000  SIM12620 

C  SIM12630 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM12640 
C  SIM12650 

7000  CONTINUE  SIM12660 

C  SIM12670 

C  SIM12680 

CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  SIM12690 


C  START  A  NEW  REGION 
C  IFCKR.EQ.IO)  GO  TO  7000 

C  KR=KR+1 

C  GO  TO  2000 


7000  CONTINUE 


C  SIM12700 
C  FIND  IF  FLOW  BECAME  COLLISIONLESS  IN  THE  WHOLE  SECTOR  SIM12710 
C  IF  POSITIVE,  STORE  F0E1,F0E2  IN  A  SEPARATE  FILE  AND  STOP  PROGRAM  SIM12720 
C  SIM12730 
C  PREPARE  DATA  FOR  THE  NEXT  SECTOR  SIM12740 
C  KR=1  SIM12750 
C  STOP  PROGRAM  IF  KS  WAS  BOUNDED  BY  THE  WALL  SIM12760 
C  KS=KS+1  SIM12770 
C  GO  TO  1000  SIM12780 
C  SIM12790 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXSIM12800 


C  SIM12810 
C  IF  THERE  IS  BACK  FLOW  (FWN1,FWN2)  NONZERO  CALCULATE  NEXT  ITERATION  SIM12820 
C  ITER=ITER+1  SIM12830 
C  KS=KR=1  SIM12840 
C  GO  TO  3000  SIM12850 
C  SIM12860 
CXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSIM12870 


WRITE(6,6099) 

6099  FORMAT(*l  DATA  FOR  TEN  MOLECULES  SPEC. 2') 

DO  3003  1=1,10 

WRITE(6,3002)  P2( 1 , I ) , P2(2, I) , P2(3, I) , P2(4, I ) , P2(5, I ) 

3002  FORMATC  •  ,4F10 . 3,  E18 . 10) 

3003  CONTINUE 
GO  TO  3009 

3004  WRITE(6,3005) 


SIM12880 

SIM12890 

SIM12900 

SIM12910 

SIM12920 

S1M12930 

3IM12940 

SIM12950 

SIM12960 
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3005 

FORMATC  NO  PLACE  FOR  ADDITIONAL  MOLECULES') 

SIM12970 

3009 

CONTINUE 

SIM12980 

STOP 

SIM12990 

END 

SIM13000 

SIM13010 

SIM13020 

SUBROUTINE  RANDU(P) 

SIM13030 

COMMON  IX 

SIM13040 

lY  »  IXX655S9 

SIM13050 

IF  CIY)  5.6,6 

SIM13060 

5 

lY  s  IY42147683647+1 

SIM13070 

6 

P  s  lY 

SIM13080 

P  =  PX.4656613E-9 

SIM13090 

IX  *  lY 

SIM13100 

RETURN 

SIM13110 

END 

SIN13120 

B.5  Program  SIMUL  -  User's  Guide 

Preparation:  Run  AXSYIl  program.  From  its  output  data  evaluate: 

ALFA  -  The  averaged  angle  for  the  continuum  breakdown 

(P  -  0.05) 

TETA  -  Flow  direction  along  the  breakdown  limit. 

Flow  parameters  along  this  line  -  Pressure,  temperature,  velocity,  mean 


free  path. 

Input  Data: 

Radius  of  nozzle  ring  R1 

Radial  size  of  a  cell  DR 

Maximum  radius  in  simulation  RP 

Angle  of  breakdown  limit  ALFA 

Flow  direction  along  (ALFA)  TETA 

Averaged  flow  velocity  along  (ALFA)  Vo 

Molecular  weight  of  each  species  Spec (1,1) 

Molecular  diameter  of  each  species  Spec  (1,2) 

Mean  free  path  along  (ALFA)  FPM 

Averaged  Temperature  (ALFA)  TEMP 

Time  Increment  DTH 

Number  of  time  Increments  NIS 

Options 

a.  Geometry 


The  program  is  designed  to  run  for  axisymmetric  ring  flow. 

For  a  two  dimensional  planar  flow  the  molecular  motion,  collisions  and 
flow  calculation  remain  unchanged.  The  flow  cross  section  remain 


unchanged.  Instead  of  the  angle  DFI  the  two  dimensional  flow  requires 
the  definition  of  the  width  of  each  region  (or  cell).  To  keep  the 
number  of  molecules  within  reasonable  computational  limits  this  size 
has  to  decrease  in  the  same  manner  as  DFI. 

The  part  of  the  program  which  has  to  be  changed  for  this  purpose  is 
lines  230  to  250. 
b.  Wall  flux  calculations: 

The  program  may  run  for  the  whole  molecular  region  resulting  the  flux 
towards  the  wall  (this  part  needs  additional  debugging).  However  in 
order  to  make  it  more  efficient  we  may  stop  the  program  at  a  sector 
where  the  flow  becomes  collisionless.  The  remaining  part  of  the  flow 
may  be  regarded  either  as  collisionless  or  if  we  define  a  much  larger 
size  of  cells  we  may  calculate  the  molecular  collisions  on  this 
bas is . 

This  part  need  additional  analysis  and  programming, 
hxecution  Commands 


Without  additional  changes  the  program  runs  under  WATFIV  compiler  using 
Che  following  command: 

WATFI  AXSYM  *  (XTYPE) 

Further  developments  will  be  required  to  run  the  program  for  each  sector 


separately  and  Che  output  intermediate  results  on  the  mass  storage. 


B.6 


The  Influence  of  the  Ambient  Gas 


The  temperature,  pressure  and  density  of  the  ambient  gas  are  shown  in 
Table  1.  Because  its  temperature  is  much  higher  than  in  the  jet  gas,  the 
thermal  velocity  of  ambient  gas  molecules  is  much  higher  than  jet  molecules, 
the  following  contribution  may  be  expected  due  to  the  ambient  gas : 

a.  Collisions  between  the  "hot**  ambient  molecules  with  the  "cold"  jet 
molecules  may  cause  an  increase  in  the  dissipation  in  the  outer 
layer  of  the  jet  and  Increase  in  the  flux  towards  the  walls. 

b.  For  higher  ambient  pressures,  those  collisions  become  rare  because 
of  the  low  number  density  therefore,  the  influence  of  the 
collisions  may  decrease. 

The  only  way  to  evaluate  the  influence  of  these  two  controversial  factors 
is  by  an  additional  simulation  program  to  be  designed  for  this  region.  The 
following  are  the  main  factors  to  be  Included  in  this  program: 

Boundary  Conditions: 

a.  The  jet  side:  FOEl  and  F0E2  obtained  from  the  last  simulated  sector 
(SIIIUL)  supply  the  number  flux,  flow  velocity  componetns  and 
temperatures.  These  parameters  are  given  for  all  points  along  the 
radius  K(  1)  at  constant  distances  DR.  In  the  low  density  domain  the 
resolution  uR  is  much  too  large  compared  with  the  expected  mean  free 
path.  A  different  mesh  has  to  be  designed  for  this  purpose.  It  is 
possible  that  one  cell  may  be  sufficient. 

b.  Far  Field  Condition:  The  boundary  conditions  where  the  gas  may  be 
regarded  "undisturbed”  by  the  jet  are  as  follows 


-  Jet  gas  molecules  are  allowed  to  go  out  the  simulation  region 
(these  molecules  will  be  regarded  as  "lost"  molecules). 

-  Ambient  gas  is  allowed  to  enter  the  simulated  region  according 
with  their  thermal  velocity  and  number  density. 

c.  Solid  Wall  Bondarles.  The  solid  wall  may  be  assumed  to  have  a 

constant  temperature  T„.  Incident  molecules  of  either  species  are 
reflected  back  from  the  wall.  Different  models  of  collisions  with 
the  wall  may  be  employed. 

-  Elastic  collisions  'specular  reglection'  (this  calculation  has  been 
Included  in  SUiUL  program) 

-  Collisions  with  ideal  heat  transfer  (diffuse  reflections) 

-  Other  models  depending  on  the  materials  and  surface  parameters. 

The  collision  with  the  wall  was  Included  as  an  internal  routine  in  the 

program  SIMUL.  If  the  general  molecular  simulation  contains  the  program 
proposed  here  the  collision  with  the  wall  will  be  omitted  from  SIIIUL  and 
become  the  core  of  the  additional  collisionless  program.  Figure  25  shows  the 
low  density  (collisionless)  region  and  its  boundaries. 


Figure  25.  The  Low  Densicy  Region  in  Che  Jet 


SUMMARY  OF  REPORT 


Algorithms  for  the  continuum  regime  and  for  the  region  where  molecular 
collisions  are  significant  have  been  developed. 

Program  AXSYM  contains  the  calculation  of  planar  Jet  flow  and 
axlsymmetrlc  ring  jet  flow.  This  program  supplies  data  for  the  limits  where 
the  continuum  approach  become  Invalid  and  molecular  approach  should  be 
employed . 

Program  SIMUL  contains  the  molecular  simulation  for  the  axlsymmetrlc  ring 
flow.  This  program  may  run  for  the  whole  molecular  region  to  result  In  the 
calculation  of  flux  towards  the  solid  wall.  For  a  more  efficient  simulation 
It  Is  proposed  to  design  an  additional  program  for  the  colllslonless  region 
where  ambient  gas  may  be  Included. 

For  the  two  dimensional  flow,  program  SIMUL  may  be  used  after  changing 
the  definition  of  the  cell  geometry. 

To  run  the  whole  program  It  may  be  required  to  make  separate  runs  for 
each  sector  and  store  results  on  the  mass  storage. 
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